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ABSTRACT 


A  computer  model  for  axially  symmetric  dielectric  radomes  in  the  near  field  of 
a  circular  aperture  has  been  developed.  This  model,  based  on  a  method  of 
moments  solution  for  bodies  of  revolution,  is  used  to  determine  the  electric  field 
pattern  of  a  linearly  polarized  circular  aperture  radiating  through  a  dielectric 
radome.  The  program  is  written  in  FORTRAN  and  can  accommodate  rotationally 
symmetric  radome  shapes.  Arbitrary  illumination  distributions  for  the  aperture  can 
be  specified,  and  phased  scanning  of  the  aperture  in  both  azimuth  and  elevation 
is  allowed.  Computational  run  time  has  been  reduced  by  taking  advantage  of 
mode  symmetry.  A  separate  program  has  been  developed  to  compute  gain,  and 
together  these  programs  are  used  to  determine  the  pattern  degradation  and  gain 
loss  due  to  the  presence  of  a  radome  in  the  near  field. 
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I.  INTRODUCTION 


A  radome  is  a  structure  designed  to  protect  an  antenna  in  its  operating 
environment.  Aircraft  radomes  must  protect  radar  antennas  from  high 
aerodynamic  forces  while  minimizing  drag,  weight,  and  interference  on  its 
electrical  operation.  The  presence  of  a  radome  can  affect  the  gain,  beamwidth, 
sidelobe  level,  and  direction  of  boresight,  as  well  as  change  the  VSWR  and 
antenna  noise  temperature.  This  thesis  involves  the  continued  development  of  a 
computer  code  that  predicts  antenna  radiation  patterns  looking  through 
aerodynamic  radomes. 

The  ability  to  accurately  predict  the  effect  of  a  radome  on  the  operation  of  an 
antenna  early  in  the  design  process  is  desirable.  Typically,  designers  of  the 
earliest  radomes  were  concerned  more  with  structural  considerations  than 
electrical  performance.  However,  after  the  near  failure  of  the  U.S.  Army’s 
"Pathfinder"  program  of  radar-equipped  B-24  "Liberator"  bombers  due  to  poor 
electrical  performance  of  their  radomes,  it  became  clear  that,"  There  was  a  serious 
need  for  developing  sound  theoretical  procedure  for  the  electrical  design  of 
radomes."  [Ref.  1]  For  most  practical  radomes  the  surface  geometry  is  too 
complex  to  yield  closed-form  mathematical  expressions  describing  their 
performance.  It  is  possible  to  predict  radome  performance  using  numerical 
solutions  of  Maxwell’s  equations.  The  program  described  in  this  thesis  is  based 
on  a  numerical  technique  called  the  method  of  moments  (MM).  The  particular 
form  of  the  solution  used  here  is  based  on  a  formulation  by  Mautz  and  Harrington 
[Ref.  7]  for  bodies  of  revolution.  Francis  [Ref.  2]  developed  a  computer  code 
LDBORMM.F  that  specialized  the  Mautz  and  Harrington  solution  to  ogive- 
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shaped  radomes.  There  are  no  restrictions  on  the  location  of  the  radome  with 
respect  to  the  antenna,  and  therefore  near-field  radomes  can  be  analyzed. 

The  goal  of  this  thesis  was  to  extend  the  capabilities  of  this  program. 
Significant  improvements  include:  taking  advantage  of  symmetry  to  reduce 
computation  time,  applying  realistic  radar  antenna  patterns  for  comparison  to 
experimental  results,  and  the  computation  of  gain  loss  due  to  the  radome. 

A.  DESCRIPTION  OF  THE  PROBLEM 

The  basic  system  analyzed  is  a  radome  consisting  of  a  body  of  revolution 
(BOR)  located  coaxial  with  a  circular  antenna,  as  in  Figure  1.1.  The  BOR 
constraint  is  satisfied  by  many  aerodynamic  radomes  and  greatly  simplifies  the 
MM  solution  ,  which  will  be  described  in  more  detail  in  following  chapters.  The 


ogive  shape  of  the  radome  in  Figure  1.1  is  the  primary  one  of  interest,  but  the 
program  can  also  model  other  bodies  of  revolution  including:  a  cone,  disk, 
hemisphere,  and  parabola,  as  shown  in  Figure  1.2.  These  were  primarily  added  for 
validation  purposes. 
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There  are  no  constraints  on  the  size,  shape  and  location  of  these  objects 
provided  they  are  coaxial  with  the  antenna.  However,  because  MM  requires  the 


Figure  1.2  Optional  Radome  Shapes 


solution  of  a  matrix  equation,  computer  run  time  and  memory  will  limit  the  size  of 
the  bodies  that  can  be  analyzed.  A  unique  feature  of  this  program  is  that  it  takes 
into  consideration  the  fact  that  the  radome  may  be  in  the  near  field  of  the  antenna 
as  in  Figure  1.3.  The  antenna  aperture  used  in  the  program  is  circular,  but  the 
program  can  be  easily  altered  to  accommodate  other  shapes.  There  are  no 
constraints  on  the  aperture  size,  scan  angle,  or  polarization. 

The  remainder  of  this  thesis  is  divided  into  three  major  sections.  Chapter  II 
covers  theoretical  development  of  the  problem  and  method  of  solution.  The  MM 
solution  is  described  in  Chapter  II  along  with  a  discussion  of  mode  symmetry  and 
how  it  can  be  used  to  advantage  in  obtaining  the  radome  current.  Chapter  II  also 
describes  how  the  total  field  is  calculated  from  the  currents.  Chapter  III  covers 
the  programs  used  for  generating  patterns  based  on  the  methods  of 
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Chapter  II.  The  first  part  describes  the  program  RADOME.F  which  calculates  the 
current  and  electric  field.  The  second  part  describes  the  program  GAIN.F  which 


Figure  1.3  Radome  Located  in  Near  Field  of  Antenna 

calculates  the  gain  of  the  system.  The  final  chapter  discusses  the  results  and 
presents  several  conclusions  based  on  the  calculated  radome  patterns. 
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II.  ANALYSIS 


A.  GENERAL  DESCRIPTION  OF  ANALYSIS 

The  objective  of  the  analysis  is  to  determine  the  radiation  pattern  of  the 

antenna  and  its  gain  in  the  presence  of  a  radome.  In  order  to  analyze  the  system, 

the  spherical  coordinates  of  Figure  2.1  are  used.  The  tote'  electric  field  in  the 
direction  of  unit  vectors  (0,o)  is  required  to  compute  the  gain  in  the  far  field.  The 

total  electric  field  is  the  sum  of  the  electric  field  due  to  the  currents  on  the 
antenna  (radiated  field),  and  the  electric  field  due  to  the  currents  on  the  radome 
(scattered  field).  The  procedure  is  to  assume  a  current  distribution  on  the 
antenna,  then  use  the  MM  to  determine  currents  induced  on  the  radome.  In 
practice,  scattering  from  the  radome  back  toward  the  antenna  will  cause  a  change 
in  the  initial  current  distribution.  This  effect  will  be  small  for  a  well-desigr.ed 
radome  and  is  ignored  in  the  analysis.  Scattering  back  toward  the  antenna  will 
primarily  affect  the  YSWR  looking  into  its  terminals. 

1.  Determination  of  Gain 

Gain  is  a  measure  of  the  ability  of  an  antenna  to  concentrate  radiated 
power  in  a  particular  direction.  The  maximum  value  of  gain  for  a  lossless  antenna 
is  the  directivity  Go,  which  is  the  ratio  of  the  maximum  radiation  intensity  to  the 
average  radiation  intensity  [Ref.  3:p  610], 


4 /r  1/(0,  <z>) 

_ i _ iman 

§U(d,<l>)dn  ' 


(2-1) 
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where  Z/(0, 0)  is  the  radiation  intensity  and  dQ  =  sin  9ddd<t>  is  the  differential 
solid  angle.  Throughout  this  paper  a  lossless  antenna  will  be  assumed  and  the 


Figure  2.1  Coordinate  System 

terms  gain  and  directivity  used  interchangeably. 

Radiation  intensity  is  the  time-averaged  power  per  unit  solid  angle,  and 
is  proportional  to  the  magnitude  of  the  electric  field  squared 

U  =  —  |e|\  (2-2) 

2n'  1 


where 
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and  Ea  is  the  electric  field  from  the  antenna,  Es  is  the  electric  field  scattered  by 
the  currents  induced  on  the  radome,  and  rj  is  the  impedance  of  free  space. 

2.  Determining  the  Electric  Field 

The  total  electric  field  is  calculated  from  the  current  distributions  on  the 
antenna  and  radome  as  illustrated  in  Figure  2.2.  A  current  distribution  Ja  on  the 
antenna  is  assumed  known.  The  currents  induced  on  the  radome  are  dependent 
on  the  current  distribution  on  the  antenna  as  shown  in  Figure  2.3,  and  this  is 
taken  into  account  in  the  method  of  moments  solution. 

Maxwell's  equations,  are  used  to  determine  the  electric  field  at  a  point  in 
space  due  to  a  current  density  J.  When  both  the  radome  and  antenna  are 
present,  the  total  current  density  is  J  =  Ja  +  j#.  A  solution  of  Maxwell's 
equations  is 

E  =  -jo)H - Lv(v.R),  (2-4) 

tope  '  ’ 

where 

fl  =  ^-JJO(r.9.0)^--— ds\  (2-5) 

and 


r 


R-R' 


(2-6) 


In  (2-4)  R  is  the  vector  potential,  R  the  position  vector  to  the  observation  point 
and  R'  the  position  vector  to  a  source  point  on  the  antenna. 


Figure  2.2  Electric  Field  from  Currents  on  Antenna  and  Radome 


Figure  2.3  Electric  Field  Impinging  on  Radome 
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B.  ANALYTICAL  TECHNIQUE 


Closed  form  mathematical  solutions  of  the  gain  and  radiation  patterns  are  only 
possible  for  a  limited  number  of  applications  that  neglect  the  presence  of  the 
radome.  In  general,  determining  the  radiation  pattern  and  gain  with  the  radome 
requires  a  numerical  solution  of  Maxwell's  equations.  The  radome  geometry 
results  in  an  unknown  current  distribution  so  (2-4)  cannot  be  solved  directly. 
However,  boundary  conditions  can  be  applied  on  the  radome  surface  which  yield 
an  integral  equation  for  the  current.  The  integral  equation  will  be  solved 
numerically  using  the  method  of  moments. 

1.  E-field  Integral  Equation  and  the  Method  of  Moments 

The  method  of  moments  (MM)  is  a  procedure  used  to  solve  an  integral 
equation  obtained  by  enforcing  (2-4)  on  the  radome  surface  [Ref.  5].  Equation  (2- 
4)  is  used  to  determine  the  electric  field  E*  incident  upon  the  surface  of  the 
radome  from  the  currents  on  the  antenna  J, .  Initially  the  surface  will  be  assumed 
a  perfect  electric  conductor  (PEC),  in  which  case  the  total  E  tangent  to  the  surface 
of  the  radome  is  equal  to  zero.  Later  a  correction  term  will  be  added  to  the  perfect 
conductor  solution  to  account  for  the  dielectric  properties  of  the  radome. 

The  boundary  condition  on  the  PEC  radome  is 

EStan  =  -Eata„.  (2-7) 


As  applied  to  (2-4), 


Mi  JJ  J%Gds  +  ^  v  v  *  JJ  j*Gds 


(2-8) 


where 
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(2-9) 


C  = 


g-Jkr 

4jtr ' 


and 


r  =  R-R 


(2-6) 


Equation  (2-8)  is  the  E-field  integral  equation  (EFIE).  The  only  unknown  is  Js, 
and  once  determined  it  is  possible  to  calculate  the  scattered  field  by  applying  (2-4) 
to  the  radome. 

To  solve  (2-8),  the  current  Js  is  expanded  into  a  series, 

3,-flA.  (2-10) 

i-i 


The  J;  are  basis  functions,  and  the  I(  are  the  expansion  coefficients  to  be 
determined. 

Mautz  and  Harrington  applied  this  technique  to  the  specific  case  of 
surfaces  S  generated  by  revolving  a  plane  curve  about  the  z-axis  [Ref. 5].  The 
plane  curve  is  approximated  by  a  series  of  straight-line  segments  connecting  Slp 
points.  The  surface  has  a  circular  cross  section  in  the  azimuth  variable  4>,  thus  the 
body  is  represented  by  a  series  of  "faceted  rings"  as  shown  in  Figure  2.4. 
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The  surface  coordinate  system  is  cylindrical  (p,  0fz),  and  1  is  a  vector 

tangent  to  the  surface  and  orthogonal  to  <p  at  every  point.  In  order  for  the  series 
expansion  of  the  current  (2-10)  to  approximate  an  arbitrary  surface  current  density 
Js,  independent  sets  of  functions  are  defined  as  follows: 


Figure  2.4  Segmented  Radome  with  Current  Basis  Functions 
Depicted 


in  t: 

in  0'. 


Tjit^J  is  a  triangle  function  extending  over  two  segments,  and  MO isa 

pulse  function  extending  over  a  single  segment.  For  an  explanation  of  these 
functions  see  [Ref.  l:pp  16-18],  Each  value  of  n  is  referred  to  as  an  azimuthal 
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mode.  The  series  expansion  of  (2-10)  in  terms  of  the  chosen  basis  functions 
becomes. 


I 


Np- 2 

I  'U 

.  J=  I 


nj 


Np-\ 

+ 


(2-12) 


In  order  for  the  series  to  be  an  exact  solution,  the  sum  must  be  extended  over 
infinite  modes  Of  course  the  series  must  be  truncated,  and  the 

maximum  number  of  modes  for  convergence  nmw  is  dependent  upon  the  antenna 
scan  angle  and  radome  size. 

When  the  series  expansion  is  applied  to  the  integral  equation  (2-8), 


(2-13) 


In  order  to  solve  for  the  unknowns  fnJ  and  lenJ,  both  sides  of  (2-13)  are  multiplied 

by  testing  functions  Ufmi  and  UJamj ,  and  integrated.  Using  Galerkin’s  Method 

[Ref.  8],  the  testing  functions  are  chosen  to  be  the  complex  conjugates  of  the  basis 
functions 


e~Jm* 

g-jm. 


(2-14) 
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Multiplying  both  sides  of  (2-13)  by  each  of  the  testing  functions  and 
integrating  results  in  a  set  of  2  Np  -  3  simultaneous  linear  equations  for  each  value 

of  H.  In  matrix  form,  these  equations  are  written  as 


(2-15) 


where  llj  is  the  excitation  vector,  Z  is  an  impedance  matrix  and 


contains  the 


unknown  current  coefficients.  Submatrices  associated  with  the  t,  and  0 
components  (2-15)  can  be  identified  as  follows: 


uy 

[*] 

In 

u • 

Zo, 

'z00' 

r 

_ 

n 

_ 

_ 

n 

n 

_ 

. 

n 

n  =  0,±1,±2 


max- 


(2-16) 


The  unknown  coefficients  are  determined  by  solving  the  matrix  equation 

M-WH  ,2-171 

The  MM  solution  for  the  unknown  current  coefficients  /!,.  and  l°„  in  vector  form 

nJ  nj 

is 


In 

[«] 

-i 

r 

Zol 

z 

u* 

_ 

n 

_ 

n 

n 

_ 

n 

_ 

(2-18) 


Once  these  are  determined,  it  is  possible  to  compute  the  electric  field  due  to  the 

current  density  distribution  Js  on  the  surface  of  the  radome  at  any  point  in  space. 
The  elements  of  each  submatrix  of  the  excitation  vector  [ll]  are 
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(2-19) 


U'  =  J-  \\uf  »Ea  1 ds 
m  n  JJ  ni 

s 

zy°  =  i  [[ur  .p  *ds, 

111  n  JJ  ni  tan  ’ 


and  the  elements  of  each  of  the  submatrices  of  the  impedance  matrix  Z  are 


=  *11 


5/  L 


S|  S/  L 

5,  5/  L 


*1 


ds'ds 

ds'ds 

ds'ds 

\ds'ds. 


(2-20) 


The  fact  that  the  testing  functions  are  orthogonal  in  0  to  due  to 
the  presence  of  the  exponentials  has  been  used  to  eliminate  elements  of  (2-20) 
with  n  *  m  [Ref.  6].  This  allows  each  mode  n  to  be  treated  independently, 
thereby  simplifying  the  matrix  inversion  problem.  Thus,  by  constraining  the 
surface  to  be  rotationally  symmetric,  the  azimuthal  dependence  of  the  current  is 
expressed  as  a  sum  over  exponentials  (a  Fourier  series)  and  each  mode  n  is 
independent. 

a.  Impedance  of  Thin- Shell  Dielectric  Radomes 

The  impedance  matrix  [z|  of  the  MM  solution  assumes  that  the 
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radome  surface  is  a  PEC,  where  (2-7)  is  applied  as  a  boundary  condition  of  (2-4). 
For  thin-shell  dielectric  radoines,  a  correction  term  is  added  to  the  impedance 
matrix  that  allows  the  modeling  of  radomes  of  arbitrary  dielectric  constant  [Ref .9], 
When  the  electromagnetic  wave  from  the  antenna  is  incident  on  a 
dielectric  radome  that  has  an  intrinsic  impedance 


n  = 


(2-21) 


different  from  that  of  free  space  [r?0j.  there  will  be  reflection  and  refraction  as 
shown  in  Figure  2.5  [Ref.  3:p.  397],  The  total  fields  are 


E  =  Ea  +  E* 

H  =  H*  +  Ha .  (2-22) 


These  fields  satisfy  Maxwell's  equations 

Vx(Ea+Es)=  -jrn0  (H8  +  H* )  (2-23) 

and, 

V  x  (H8  +  H* )  =  jcoe0( V  +  Es )  +  jco(e  -  )(E8  +  Es ) ,  (2-24) 

where  the  second  term  exists  only  when  there  is  material  present  with  a 
permittivity  different  than  free  space.  This  term  is  a  polarization  current  J,  where 

J  s  jco[e  -  e0 )E.  (2-25) 
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Maxwell's  equations  are  also  satisfied  by  the  field  from  the  antenna  (Ea,Ha). 
Since  (E.flj  and  (Eaf  Ha)  satisfy  Maxwell's  equations,  by  superposition  and  (2-22) 
(Ef  .M*  j  also  satisfy  Maxwell's  equations.  Therefore 


V  x  H*  =  jcoV  +  J. 


(2-26) 


Figure  2.5  Discontinuity  of  Permittivity  for  Dielectric 
Radomes 


For  a  radome  consisting  of  a  thin  dielectric  shell,  the  polarization 
current  has  a  large  tangential  component  and  small  normal  component.  If  only  the 
tangential  component  is  significant,  then  (2-8)  can  be  modified  to  include 
dielectric  properties  of  the  radome  as  follows: 
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(2-27) 


§  a  _ 

Man  — 


J0)(i 


=5 — -  +  jeon  f  f  J%Gds  +  -J-  V  V  .  ff  jsGds 
-£0)t  V  7 


The  first  term  on  the  right  side  is  the  impedance  loading  term;  the  integro- 

differential  portion  remains  unchanged  from  the  PEC  assumption. 

After  multiplying  both  sides  of  (2-27)  by  the  testing  functions 
[Wlm , UJ^i )  and  integrating,  the  impedance  matrix  jzj  of  (2-15)  becomes 


[ZHZmmMZJ.  (2-28) 

where  the  elements  of  [ZMM  j  are  given  by  (2-20),  and 

[*".], -P* 

S 

[^i»]  =  JJ Hfj  •  J °%jnZLdS  -  0,  (by  orthogonality  of  t,  and  o) 

[Z^j  =  jj &,*  •  Jf%jnZLdS  =  0,  (by  orthogonality  of  o.  and  t) 

S 

[z?.}ij=llUV'S,J.ZldS.  (2-29) 

5 


and 


1 

Ja)(e-eB)t 


(2-30) 


The  thickness  t  in  the  impedance  loading  term  is  necessary  to  convert  the  volume 
current  density  j,  of  Figure  2.5  to  a  surface  current  density  Js. 
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Solving  (2-17)  yields 


I 


(2-31) 


The  radome  impedance  ZL  can  be  rewritten  in  terms  of  dielectric  constant  er , 
thickness  to  wavelength  ratio  nr  .  and  the  loss  tangent  tan  8  [Ref.  2:pp  45-46], 


ZL 


60 


n. 


j(ef  —  l)  +  £r  tan S 


(2-32) 


Note  that  both  (2-29)  and  (2-32)  allow  the  modeling  of  radomes  whose 
permittivity  may  be  complex,  e  =  e’  -  je"  . 
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2.  Mode  Symmetry 

When  the  surface  current  density  Js  is  expanded  into  a  series  as  in  (2- 
12),  the  sum  is  extended  over  azimuthal  modes  n  =  0,±1,±2,  •  ,±nmaM.  The 
mode  index  appears  in  the  evponential  term  of  basis  functions  (2-11)  and 
testing  functions  (2-14).  The  exponential  factor  can  be  written  as, 

ejn  =cosino)  +  jsin[n0).  (2-33) 

For  each  n  there  is  a  current  vector  ln  given  by  (2-18).  Mode  symmetry  refers  to 
a  relationship  between  the  positive  and  negative  inode  vectors  l„  and  l_„.  If  mode 
symmetry  exists,  then  only  one  matrix  equation  needs  to  be  evaluated  rather  than 
two.  This  greatly  reduces  the  computational  effort. 

Recall  the  definition  of  the  excitation  vector  [Oj 

K  =  -  jj  ufm,  •  i'Jds 

s 


The  electric  field  from  the  antenna  is 

Ea  =F*(Hf0,o)R  +  £‘*(/?f0,0)0  +  £7(»,0f0)d,  (2-34) 

which  on  the  radome  surface  can  be  expressed  in  terms  of  the  two  orthogonal 
components  f  and  <2> 

E*  =£*(<,«>)»  +c:(t,0)i.  (2-35) 
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In  the  far  field  of  the  antenna,  the  R  component  in  (2-34)  decays  to  zero  and  Ef 
consists  only  of  £*.  Therefore  as  shown  in  Figure  2.6, 


E‘{t,o)  =  E‘,(t)sin$.  (2-36) 


X 


Figure  2.6  Components  of  the  Antenna  Electric-Field  on  the  Radome 

Plugging  these  far-field  components  into  (2-19) 

Wm  =  J  T^Ef(ty(t)dt  J  cosp[cos(/70)  +  jsin(n0)]fJ<t> 

f*l  *  0a0 

=  j  ^Ea9(t)i)(t)dt  J  sin0[cos(n<>)  + Jsin(/?0)]d0.  (2-37) 

'  o>a  '  j 
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Only  the  0  integration  is  dependent  on  n.  Expanding  the  integral  in  <p  and 
explicitly  including  the  sign  oi  n  gives 

2  it 

J  coso  cos(±no)  +  7'sin(±n0)]£/0  = 

•=0 

2k  2k 

j  cosocos(rn0)£/d  +  j  Jcos0sin(±o<?)rfp,  (2-38) 

#=a  «=8 

where 

2- 

jcososin(±/?0)£/0  =  0,  (2-39) 

<5  =  0 

due  to  the  fact  that  the  integrand  is  an  odd  function  of  <p ,  and  the  integration  is 
over  one  period  0  <  0  <  2n.  By  trigonometric  identity 

2 *  2  it  2it 

J  cos0cos(±n0)t/0  =  J  cosm ±n - +  J  cosl(±/?  +  l)p  da>.  (2-40) 

«*0  =8  «=0 

The  right  hand  side  is  only  non  zero  if  H  =  ±1.  Furthermore,  because  cosine  is  an 
even  function,  the  term  lfni  is  independent  of  the  sign  of  n.  Thus 

(2-41) 

For  the  second  equation  of  (2-  *7 ),  the  sign  of  Z/®(  is  dependent  on  the  integral 
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(2-42) 


J  sin 0  cos(±#i0)  +  jsin(±/ty)]tf0  = 

#=8 

2k  2k 

J  sin  o  cos(±n0}d0  +  j  Jsin0sin(±/70)</0, 


0*0 


where 


2.- 

J  sin0cos(±/70)tf0  =  0, 

o=0 


(2-43) 


by  symmetry.  By  trigonometric  identity 


2 k  2k  2k 

J  sin0sin(±/?0)</0  =  |  cos f(±n- 1)000  -  J  cos  (±/i  +  i)0|f/<?.  (2-44) 

«=fl  o=8  »=8 


As  can  be  seen  in  (2-44),  the  excitation  vector  U*nj  is  dependent  on  the  sign  of  n 
due  to  the  presence  of  the  negative  sign.  Therefore 

U°  =  -U°.  (2-45) 

-nj  +nj  v  ' 


In  the  near-field  of  the  antenna  the  0  dependence  is  more  complicated 

than  the  sine  and  cosine  functions  in  (2-36).  However,  it  has  been  verified 
numerically  that  F*  and  FJ  ate  always  even  functions  of  0  and  F*  is  always  odd. 

This  conclusion  assumes  that  a  symmetric  amplitude  and  linear  phase  distribution 
excites  the  antenna. 

In  light  of  the  above  results,  the  elements  of  the  impedance  matrix  [ZmM] 

of  (2-18)  have  the  following  relationship  between  positive  and  negative  modes: 
(Ref.  6] 
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7n 

u 

7 i  * 

U 

\-T* 

i; 

l*1 

*-MM  (-») 

ij 

-JO  0 

ij  . 

\l*° 

U  . 

.  (2-46) 


By  inspection  of  (2-30),  the  elements  of  the  load  impedance  matrix  |Ztj  are 
independent  of  the  mode  /?,  because  UJjZ  and  J'®,  are  complex  conjugates  of  each 
other.  Therefore, 


\j" 

\e\ 

iv;  .1 

w  ' 

[Zt]  = 

L  L(-")J 

u  L  J 

L  L(+")Jy 

\a] 

h-li 

\b] 

|Z"  ,1 

i  \ 

[  \ 

L  J 

(2-47) 


Rewriting  (2-3 1 )  results  in 

l'WvM 

where 

[V]  =  [7M„+Zl]"'.  (2-49) 

By  taking  the  inverse  through  partitioning  [Ref.  8]  it  can  be  shown  that  the  mode 
symmetry  survives  the  inversion,  and 


V" 

T~n(j . 

[y,# 

-nij 

V" 

■  *mj 

-Y" 

*mj 

V*' 

-nij 

uot 

-nij 

-V  •' 

+nij 

VM 

L 

(2-50) 


By  multiplying  out  (2-48), 


23 


(2-51) 


Equations  (2-41),  (2-45),  and  (2-50)  relate  the  positive  and  negative  modes  of  (2- 
51),  as  follows: 


t.  -vvc-  »^(-":w) 

(2-52) 


This  is  an  extremely  important  result.  Using  mode  symmetry  reduces  the  number 
of  calculations  required  to  find  the  current  coefficients  l£*  by  nearly  half  in  most 

cases. 
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3.  Measurement  Matrices  for  the  Scattered  Field 

Once  the  current  *js  is  obtained  the  scattered  field  E*  can  be  determined 
by  summing  the  contributions  of  the  individual  current  elements  in  the  far-field  as 
shown  in  Figure  2.7.  This  is  accomplished  by  use  of  a  measurement  vector  with 
elements 

Rm  =  Jj  (R  +  0  +  0)jB'Jkr  •  Jst/S.  (2-53) 


E* 

In 

far-field 


Figure  2.7  Scattered  Field 

In  (2-53),  a  unit  radiated  plane  wave  is  weighted  by  the  current  element  at  each 
point  on  the  surface.  Neglecting  the  R  component  in  the  far-field,  the 
measurement  vector  separated  into  6  and  0  components  is 
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(t J.)*'.  (2-54) 

s 

where  (/,£/,  and  IV  are  the  direction  cosines, 

u  =  sindcosd 
u  =  sintfsintf 

iv  =  cos  0.  (2-55) 


Thus,  the  components  of  the  scattered  field  at  a  far-field  point  are  the 
superposition  of  all  surface  segments 


(71 


(2-56) 


4.  Closed-Form  Far-Field  Antenna  Pattern 

When  computing  the  antenna  contribution  to  the  total  field  in  the  far 
zone  the  higher  order  £  terms  can  be  neglected,  and  the  radiation  integral  can  be 

reduced  to  a  closed  form.  The  source  of  this  field  is  assumed  to  be  a  uniformly- 
excited  circular-aperture  antenna  scanned  to  a  direction  (#S,0S).  Uniformly 

excited  refers  to  a  constant  amplitude  and  a  linear  phase  distribution  of  current  on 
the  antenna.  The  R ,  9  and  0  components  in  the  far-field  are 
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(2-57) 


where  Lc ,  L6 


»  _  -jk  ** 
B  4*  r 

ca  _  -Jk  e-Jkr 

0  4/r  r 

and  are  given  by 


(£,  +  „*,) 
(L,  +  i 1*.). 


L  =  JJ Me~Jkr  cosv  ds  =  0,  (M  =  0X 

S' 

N  =  JJ  jeJkrco%¥ ds'  =  JJ  (hJ„  +y  Jy  +  zJz)e-Jkr'C0S't'ds'  (2-58) 

S  S' 


where  r#  and  are  as  shown  in  Figure  2.8  [Ref.  10:pp.  455].  If  the  antenna  is 
x-polarized,  then  Ny  =  Nz  =  0  and 


N„  =  JJ^e-^^i/s'. 

s' 


(2-59) 


The  equivalent  in  spherical  coordinates  is 

Ne  =  Jj  JH  cos  e cos  0e~Jkr'cosvds' 

s 

Nt  =  JJ  JH  sin  0eJkr  coS¥  ds' .  (2-60) 

s 
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The  phase  of  the  source  point  relative  to  the  field  point,  r'CQS  yr\  is  a  length  that 
can  be  written  in  Cartesian  coordinates  as 

r'cos  y/'  =  H'u  +  y'u  +  z'uj,  (2-61) 


Figure  2.8  Circular  Aperture  Antenna  Source  Point 


where  u,  u,  and  IV  are  the  direction  cosines  (2-55).  Including  beam  scanning 
adds  an  exponential  phase  factor  to  the  current 


j  _  j  eJk\nv,*yv,*zw,) 


(2-62) 
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where  us ,  us  and  UJS  are  the  direction  cosines  of  the  scan  angles  (#*.0,)  which 

also  obey  (2-55).  Jg  is  the  magnitude  distribution,  assumed  to  be  constant,  and 
(tf'.t/'.Z')  are  the  coordinates  of  a  source  point.  For  an  antenna  lying  in  the  x-y 

plane,  Z'  =  0 .  Therefore,  N9  can  be  rewritten  as 

Ne  =  J0  cos  # cos  o\\ eJk{HV,+y u,]e~Jk{H u*yu]ds' . 

=  COSHCosajj  (2-63) 

s 

By  setting  H'  =  p'COS0'  and  y'  =  p'  Sinp',  the  exponent  term  of  (2-63)  becomes 

kp '  cos  c >'(ut  -  u)  +  sin  p'(us  -  u)  .  (2-64) 

Defining  the  variables 

fl  =  (u,-u) 

B  =  (us-u),  (2-65) 


and  considering  these  to  be  orthogonal  vectors  as  in  Figure  2.9,  it  is  possible  to 
rewrite  (2-64)  as 
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From  Figure  2.9  the  factors  containing  R  and  B  in  the  brackets  can  be  replaced 
with  COS  /  and  Sin/ 


kp'^R2  +  B2  cos  </>' cosy  +  sin0'sin/]. 


(2-67) 


By  trigonometric  identity. 


kp'^R2  +  B2  cos(tf>'  -  /)] . 


(2-68) 


Ir2  +  b 


Figure  2.9  Representation  of  R  and  B  as  vectors 


Using  this  form  of  the  exponent  in  the  integral  of  (2-63),  gives 


Ng  =  J0  cos^  cos  o  J  JeJ*pVfl2+®ico^*'Vtfp'd0'’ 


o'»0  p'=0 


(2-69) 
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which  conforms  to  the  definition  of  a  Bessel  function.  The  same  integration  is 
performed  in  the  second  equation  of  (2-60),  and  the  same  Bessel  function  results. 
Finally, 


Ng  -J0  cos 6) COSO  na2  2 


J,[ka^jR2  +  B2 
ka^lR2  +  B2 


N,  =  -J0  sino  na 2  2 


J,[  ka^JR2  +  B2 


ka^lR2  +B2 


(2-70) 


where  J,  is  a  Bessel  function  of  the  first  order.  The  far-field  equations  (2-57) 
become 


£e  =  -co$0co$<p  na2 

4  n  r 


l^kailR 


2  ,  n2 


ka^R2  +  B2 


E0  =  J„^-- — sino  na2 
4  n  r 


2J,(  fcaVfl2  +  B2 


kaiR2  +  B 2 


(2-71) 


where 


/?  =  us  -u  =  sine,  cos0s  -sinacos^ 
B  =  us-u  =  sine,  sin<z>,  -  sin  6  sino. 


(2-72) 
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V/?2  +  02  =  ^sin2  0  +  sin2  ds  -2sindssin0cos(0-0s).  (2-73) 


Using  equations  (2-71)  through  (2-73)  are  computationally  less 
intensive  than  a  numerical  integration.  These  equations  were  employed  for  both 
the  far  field  pattern  and  gain  computations. 

5.  Near  Field  of  the  Antenna 

In  the  previous  section,  the  far-field  of  a  circular  aperture  antenna  was 
determined.  When  a  radome  is  located  in  the  near-field,  the  integrals  cannot  be 
reduced  to  closed  form  and  they  must  be  evaluated  numerically.  The  near  field  is 
generally  considered  to  be  the  region  surrounding  the  antenna  where  r  <<  A ,  or 
kr  «  1  [Ref.  10:pp.  106]. 

The  complete  expressions  for  the  near-field  of  a  circular  aperture 
antenna  excited  by  a  current  polarized  in  the  x-direction  are 


=  JJ  [M.  +  (*  ~  YJkrdH'dy' 

=  JJ  [(*  -  -  IT  )M.  YJkrdH'dy' 

s' 

=  &  If  [l  ■ «  -  »')(*  -  Z')c;  J.  y*'dK  dy . 


(2-74) 


where 


k2r2  - 1  -  jkr 
r5 

3  +  5  jkr  -  k2r2 


(2-75) 
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and 

r  =  -\j(n  -x'f +(y- yf  +  (z  -  rf .  (2-76) 

The  primed  coordinates  designate  source  points  and  the  unprimed  an  observation 

point  on  the  radome  as  in  Figure  2.10.  [Ref.  11]  S'  is  the  area  of  the  antenna 
aperture.  Since  the  radome  is  rotationally  symmetric,  the  components  lEMEyEz) 

can  be  described  in  terms  of  cylindrical  coordinates  using  the  transformations 


z 


Figure  2.10  Rectangular  Coordinates  for  Circular  Aperture  Antenna 

H'  =  p'COSp' 

y'  =  p'sinp' 

dn'dy'  =  p'dp'dp.  (2-77) 

For  an  antenna  scanned  to  (flj.P,)  the  current  distribution  is 
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Since  the  expression  for  the  excitation  vector  requires  spherical  components. 
[EHEyEg  j  are  converted  using  a  coordinate  transformation 
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E„  =  Eh  s\i\6coso  +  Ey  sindsin $  +  Ez  cose 
E =  EH  cos  6COS0  +  Ey  cos0sin0-£zsin0 
Eg  =-EHs\T\0  +  Eycos0. 


(2-81) 


These  are  the  equations  programmed  to  compute  the  E-field  impinging 
on  the  radome.  The  magnitude  of  the  current  distribution  is  allowed  to  vary  as  a 
function  of  p' .  Various  antenna  sidelobe  levels  can  be  modeled  by  simply 
changing  the  function  J0(p') 
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III.  PROGRAMS 


The  mathematical  formulas  developed  in  the  previous  chapters  were 
computer  programmed  using  FORTRAN.  The  pattern  and  gain  calculations  are 
done  in  two  separate  programs:  RADOME.F  and  GAIN.F ,  both  of  which  run  on  a 
Sun  SPARCstation™  under  UNIX.  A  brief  description  of  each  follows. 

A .  DESCRIPTION  OF  THE  PROGRAM  RADOME.F 

This  program,  developed  by  Francis  [Ref.  2],  determines  the  series  expansion 
coefficients  of  the  current  «JS  induced  on  the  surface  of  the  radome.  It  also 
calculates  the  total  electric  field  surrounding  the  system,  and  generates  radiation 
patterns  using  MATLAB™.  Originally  called  LDBOR.F,  (an  abbreviation  of 
"Loaded  Bodies  of  Revolution")  in  Ref.  2,  it  has  been  modified  to  take  advantage 
of  symmetry,  and  also  allows  the  antenna  mainbeam  to  scan  in  both  0  and  0  . 

Figure  3.1  is  a  flow  chart  of  the  program.  To  run  the  program,  data  is  entered 
interactively,  except  for  a  data  file  "gaus/gaus(ITH)"  containing  the  integration 
constants  for  Gauss  Quadrature,  which  is  read  in  the  program.  All  integrals  are 
evaluated  using  this  technique  (see  Appendix  A.).  Table  I  is  a  sample  input. 

After  all  data  is  entered,  the  coordinate  points  of  the  surface  are  computed 
and  stored  in  the  arrays  RH(  =  p )  and  ZH(  =  Z ).  The  surface  impedance  of  each 
segment  is  stored  in  ZLO.  The  elements  of  the  MM  impedance  matrix  are 

computed  in  subroutine  ZMAT,  and  the  excitation  vector  elements  in  GENEX. 
These  provide  [ZMMj  and  [uj  in  equation  (2-17).  Subroutines  DECOMP  and 

SOLVE  solve  the  matrix  equation  usine  Gaussian  elimination.  Subroutine 
PLANE  calculates  the  measurement  matrices,  which  are  subsequently  used  to 
compute  the  scattered  field  from  the  radome. 
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Two  antenna  radiation  subroutines  are  included.  CIRCRTP  calculates  the 
electric  field  components  and  Ee  according  to  equations  (2-81).  It  is  used 

by  subroutine  GENEX  to  obtain  the  excitation  vector.  The  second  antenna 
subroutine  ANTFF  assumes  that  the  observation  point  is  in  the  far-field.  This  is 
used  in  the  main  program  to  compute  the  antenna  contribution  to  the  radiation 
pattern  E*  in  (2-3). 

Output  includes  an  ASCII  file  "outldbor"  which  lists  all  the  input  parameters, 
radome  geometry  and  impedance,  and  the  pattern  magnitude.  The  current 
coefficients  are  written  on  a  file  "curcoefsdat"  so  that  further  patterns  can  be 
generated  without  recomputing  the  coefficients.  Finally,  MATLAB™  formatted 
files  are  generated  for  plotting  the  d  and  0  components  of  the  electric  field. 
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TABLE  I.  SAMPLE  INPUT  RADOME.F 


-,'v-asa 


ENTER  A  LETTER  TO  INDIC  ATE  WHICH  ITERATION  THIS  IS  A 
SELECT  BOR  GEOMETRY  BY  NUMBER 


NUMBER 


BOR 

OGIVE 

SPHERE 

(ONE 

DISK 

PARABOLA 


ENTER  SURFACE  CURVATURE  (wavelengths)  7. 

ENTER  ZPRIME,WHERE  CURVATURE  STARTS  (wavelengths)  0 
ENTER  BASE  RADIUS  (wavelengths)  3 

ENTER  FILENAMES  gaus### 


ENTER  THE  FILENAME  IN  T  ( NT) 

ENTER  THE  FILENAME  IN  PHI  (NPHI) 

ENTER  THE  FILENAf  IE  IN  X  AND  Y 

ENTER  THE  PLOTTING  INCREMENT  IN  DEGREES 

ENTER  THE  HIGHEST  MODE 

ENTER  PHI  (observation)  IN  DEGREES 

ENTER  SCAN  ANGLES  IN  DEGREES:  THETA,  PHI 

ENTER  COMPLEX  IMPEDANCE.  OHMS:  (Real,  Imag) 

ENTER  THE  ANTENNA  RADIUS  (wavelengths) 


gaus2 

gaus48 

gaus20 


30,90 

(0,-1700) 
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/ Enter  Geometry  Selections?  / 


NP=  #  points  on  curve  generating  BOR 
ZH.RH=C oordinates  of  BOR 
DM=Diameter  of  Parabolic  Radome 
FOD=F/D 

ZP=Location  of  antenna  z' 

B  ASE*Base  diameter  of  antenna 


Enter  filenames  for  Gauss  in  t 


■ . .  tfiii  iiiiiia  i 


Figure  3.1  RADOME.F  Flow  Diagram 
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Read  ITH.SELECTION.B  ASE 
NP.MT,  MP.N.NROWS , 
MODES,NBLOCK,THS.I’HS. 
ARAD.PHI.IMP.ZP.RS.R  H. 
ZH.ZLO  from  "curcoefsdaK  ITH)" 


Calculate  IT.MHI.MP.MT.N 


3 


—  — ^  Prepares  for  E-field  calculation 


Figure  3.1  RADOME.F  Flow  Diagram  (Continued) 


3 


Exit - - 1  v  Loop 

»— (  DO  52  I=l,NPrT1v 


C^ICALC=<p 

Call  ZLOAD 
returns  ZL 


Subroutine  ZLOAD  calculates 
elements  of  [Zl] 
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m 


m 
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Figure  3.1  RADOME.F  Flow  Diagram  (  Continued) 
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Figure  3.1  RADOME.F  Flow  Diagram  (  Continued) 


B.  GAIN.F 


This  program  determines  the  gain  of  the  system  by  solving  (2-1).  In  order  to 
perform  the  integration,  the  current  coefficients  from  RADOME.F  must  be 
provided  so  that  the  radome  scattered  field  can  be  computed.  Input  consists  of 
specifying  an  iteration  letter,  and  selecting  the  number  of  divisions  in  the  Q  and  0 
integration.  Input  data  is  then  read  from  a  file  "curcoefsdat"  that  is  appended 
with  the  iteration  letter  chosen. 

Figure  3.2  is  a  flow  chart  of  the  program,  and  Table  II  is  a  sample  input. 
GAIN.F  contains  several  of  the  same  subroutines  that  occur  in  RADOME.F.  They 
include  PLANE,  CIRCRTP,  and  ANTFF. 

TABLE  II.  SAMPLE  INPUT  GAIN.F 


ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS  A 
ENTER  NUMBER  OF  DIVISIONS  IN  THETA  DIRECTION, 

NDIVT  =  ?  2 

ENTER  NUMBER  OF  DIVISIONS  IN  PHI  DIRECTION,  NDIVP  =?  3 
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o 


Figure  3.2  GAIN.F  Flow  Diagram  (  Continued) 
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X' 


Write  gain  values  to  "gainoutt  ITH).m 


X 


"gainout(lTH).tn"  is  formatted  for 
input  to  MATLABTM 
for  plotting  purposes 


Figure  3.2  GA/N.F  Flow  Diagram  (  Continued) 
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IV.  RESULTS  AND  CONCLUSIONS 


A.  VALIDATION  AND  CALCULATED  RESULTS 

Part  of  the  continued  development  of  this  research  is  validation  of  the 
programs.  Test  cases  were  run  to  see  how  the  programs  performed  against 
known  results.  Test  radomes  ranged  from  nearly  nonexistent  disks  located  at 
large  distances,  to  perfectly  conducting  parabolas  located  in  the  near-field.  E- 
plane  and  H-plane  patterns  defined  in  Figure  4.1  were  plotted  to  determine  the 


E-plane  H-plane 


Figure  4.1  Orientation  of  E-plane  and  H-plane 


effects  of  the  radome.  The  antenna  pattern  without  the  radome  was  plotted 
(dashed  line)  along  with  the  pattern  of  the  complete  system  (solid  line)  to 
illustrate  the  effects  of  radome  scattering. 
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1.  Code  Validation:  Scattering  from  a  Small  Disk 

The  pattern  of  an  antenna  radiating  in  the  presence  of  a  negligibly-small 
circular-disk  located  at  a  large  distance  was  computed.  Since  the  radome 
scattered  field  is  negligible,  the  predicted  output  is  that  of  a  uniformly-excited 
circular-aperture  antenna.  Figure  4.2  shows  the  plots  of  the  E-plane  pattern  for 
three  different  radius  antennas  (.5  X.  1.0  X  and  2.0  X).  The  plots  appear  as  single 
lines  with  no  scattering.  The  gain  as  computed  by  GAIN.F.  was  compared  to  the 
theoretical  gain  of  a  circular  aperture  with  a  uniform  current  distribution 


where  a  is  the  aperture  radius  [Ref.  10:pp.  483 1.  Theoretical  gains  for  the  .5  X. 
1.0  X  and  2.0  X  radius  apertures  are  9.94.  15.96.  and  21.98  dB  respectively.  In 
each  case  the  numerical  gain  is  within  1.5  dB  of  the  theoretical  gain. 

2.  Code  Validation:  Perfectly  Conducting  Paraboloid 

Programs  have  been  written  that  use  the  MM  solution  for  bodies  of 
revolution  developed  by  Mautz  and  Harrington  to  predict  the  patterns  and  gains 
of  reflector  antennas  [Ref.  7|.  Assuming  that  the  radome  is  a  perfectly 
conducting  parabola,  the  output  of  RADOME. F  and  GAIN.F  can  be  compared 
with  the  output  of  these  existing  programs. 

Figure  4.3  is  a  test  case  where  the  radome  is  a  parabola  with  surface 
impedance  Zt  -  0+  J0  Q.  The  dimensions  of  the  radome  matched  those  of  a 
reflector  which  was  analyzed  by  a  separate  computer  program  that  calculates 
both  pattern,  and  gain.  The  pattern  and  gain  for  the  perfectly  conducting  radome 
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computed  by  RADOME.F  and  GAIN.F  are  identical  to  that  computed  for  the 
reflector. 


Theta  Angle  (deg) 

Figure  4.3  Perfectly  Conducting  Paraboloid 
3.  Dielectric  Radomes 

To  determine  the  effect  of  various  radome  materials  on  antenna 
performance,  radomes  with  several  different  surface  impedances  were  analyzed. 
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Surface  impedances  were  calculated  using  (2 o2)  for  a  range  of  er  and  tan  S 
typical  of  radome  materials.  The  real  and  imaginary  components  of  ZL  as  a 
function  of  er  and  tan  8  are  plotted  in  Figure  4.4.  As  shown  in  the  figure,  the 
reactance  is  essentially  independent  of  tan  <5  over  the  range  0  <  tan  S  <  0.1. 
and  the  real  component  of  ZL  is  strongly  dependent  on  tan  5,  particularly  in  the 
region  of  low  er . 

A  test  radome  consisting  of  an  ogive  with  dimensions  given  in  Table  I 
and  parameters  nx  =  l/28  and  tan  <5  =  0.01  was  analyzed.’  The  patterns  were 
calculated  for  dielectric  constants  ranging  2  <  er  <10,  and  the  results  are 
summarized  in  Table  III.  The  E  and  H-plane  patterns  of  an  ogive  radome  with 
ZL  =34-  J1700  Q  are  shown  in  Figures  4.5  and  4.6.  The  impedance 
corresponds  to  a  radome  made  out  of  a  material  of  er  =2.0  and  tan  <5  =  0.01. 
The  scattering  is  symmetric  with  respect  to  0  for  the  unscanned  antenna,  and  the 

backlobe  level  is  -45  dB  below  the  peak.  The  scanned  case  has  higher  sidelobes. 
as  expected,  with  a  lobe  of  -30  dB  at  an  angle  270°.  The  gain  loss  due  to  the 
radome  for  the  unscanned  case  is  approximately  0.1  dB,  and  for  the  scanned  case 
0.2  dB.  The  discontinuity  in  the  H-plane  patterns  at  0  =  90°  is  due  to  the 
antenna  model.  When  0  =  90°  the  field  is  <J>  polarized,  and  according  to  (2-71), 
£*  has  no  COS  8  factor  to  force  the  pattern  to  0  at  0  =  90°.  Note  that  the 

radome  ohmic  loss  has  not  been  included. 

Increasing  the  dielectric  constant  results  in  a  constant-shape  scattering 
pattern  with  increasing  magnitude.  Figure  4.7  is  the  E-field  pattern  for  er  =10. 
which  corresponds  to  some  ceramics.  It  is  apparent  that  the  scattering  by  the 


*  These  dimensions  are  not  intended  to  match  any  particular  existing 
radome. 
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radome  in  the  rear  hemisphere  is  significant,  and  the  gain  loss  is  nearly  4  dB. 
Table  III  shows  the  gains  with  and  without  the  radome  for  2  <  er  <10.  Figure 


_ Loss  tangent  =  0.000  - 

+++  Loss  tangent  =  0.00 1 
—  Loss  tangent  =  0.010 
Loss  tangent  =  0. 100  - 


± . fc  •:  ~~  +  ~  ~  — 1»  -  -  ■» - -r  -  - + 

2  4  6  8  10 

Relative  Dielectric 
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Figure  4.4  Surface  Impedance  ZL  as  a  Function  of  er  and  tan  5 
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4.8  is  a  plot  of  the  gain  loss  relative  to  no  radome  as  a  function  of  er.  For 
reference  the  Fresnel  reflection  coefficient  for  a  planar  interface  is  also  shown. 
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TABLE  in.  GAIN  LOSS  FOR  AN  OGIVE  RADOME 


£r 

(°) 

=  ) 

[tan<5  =  0.01  j 

Gain  (dB) 
WithAVithout 
Radome 

9,  =0°  (t>,  =0°) 

Gain  (dB) 
With/Without 
Radome 

ds  =  30°  (a,  =  o°) 

2.0 

0-jl700 

22.21/22.32 

21.61/21.78 

2.0 

34-jl700 

22.21/22.32 

21.61/21.78 

4.0 

7.5-j560 

21.42/22.32 

20.90/21.78 

5.0 

5.3-j420 

20.91/22.32 

20.46/21.78 

7.0 

3.3-j280 

19.95/22.32 

19.47/21.78 

8.0 

2.7-j240 

19.51/22.32 

18.96/21.78 

9.0 

2.4-j210 

19.04/22.32 

18.47/21.78 

10.0 

2.1-jl90 

18.59/22.32 

18.00/21.78 
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Figure  4.5  E-Plane  Pattern  of  Antenna  with  an  Ogive  Radome 

Zt  -  34  -  j\  700  Q 
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Theta  Angle  (deg) 


Theta  Angle  (deg) 

Figure  4.6  H-Plane  Pattern  of  Antenna  with  an  Ogive  Radome 

Zt  -34- J1700  Q 
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Figure  4.7  E-Plane  Pattern  of  Antenna  with  an  Ogive  Radome 

Zt  -  2.1  -  j\  87  Q 
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The  gain  loss  appears  to  be  linear  with  er .  This  is  probably  because  the  base  of 
the  ogive  is  open,  which  allows  reflected  waves  to  exit  the  radome.  If  a  back 
surface  was  present,  these  waves  would  be  reflected  again  causing  additional 
gain  loss. 

a.  Conical  Radome 

To  determine  the  effect  that  radome  shape  has  on  scattering,  a 
conical  radome  was  analyzed.  The  length,  width,  and  impedance  of  the  cone  was 


2  4  6  8  10 

Relative  Dielectric  of  Radome 

Figure  4.8  Gainloss  Relative  for  an  Ogive  Radome  as  a  Function  of  er  (E- 

Plane  Scan) 

the  same  as  the  ogive  radome  discussed  previously.  A  conical  radome  presents  a 
relatively  flat  surface  that  could  possibly  cause  more  focused  scattering  than  a 
curved  shape  like  an  ogive.  Ray  tracing  of  Figure  4.9  shows  that  the  reflections 
from  the  radome  will  add  coherently  for  the  cone,  and  non-coherently  for  the 
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ogive.  However,  this  assumes  that  the  radome  is  in  the  far  field  of  the  antenna. 
The  E-plane  pattern  is  shown  in  Figure  4.10. 


Figure  4.9  Scattering  from  Conical  Radome  vs.  Ogive 

A  comparison  of  Figures  4.5  and  4.10  demonstrates  that  there  is 
even  less  scattering  from  the  conical  radome  in  the  rear  hemisphere.  For  the 
scanned  case,  the  conical  radome  has  a  distinct  lobe  at  d  =  280°  due  to  the 
specular  reflection  of  the  main  lobe  from  one  side  of  the  radome.  This  is  referred 
to  as  a  reflection  lobe. 


Figure  4.10  E-Pane  Pattern  of  Antenna  with  a  Conical 
RadomeZt  -  34  -  jl  70B  Q 
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B.  CONCLUSIONS 

This  research  has  resulted  in  two  computer  codes  that  can  be  used  to  predict 
the  pattern  degradation  and  gain  loss  due  to  an  axially  symmetric  radome  in  the 
near  field  of  a  linearly  polarized  circular  aperture.  The  computed  patterns  show 
that  the  gain  loss  is  slightly  greater  than  the  reflection  loss  for  a  planar  interface 
with  the  same  dielectric  constant  as  the  radome.  This  is  approximately  true  for 
low-loss  radomes.  For  lossy  radomes  the  impedance  has  a  very  strong  affect  on 
scattering,  while  radome  shape  has  less  effect.  Scanning  the  antenna  results  in  an 
asymmetric  scattering  pattern  and  generally  higher  sidelobe  levels. 

The  computed  gain  compares  favorably  with  values  computed  by  closed  form 
solutions  derived  from  theory,  and  the  gain  of  the  system  has  the  expected 
decrease  as  the  level  of  scattering  increases.  For  the  specific  case  of  the  ogive 
radome,  there  is  a  linear  dependence  of  the  gain  loss  to  the  relative  dielectric  of 
the  radome.  Note  that  the  gain  calculations  presented  here  do  not  include  ohmic 
loss  inside  of  the  radome. 

Improvements  in  the  original  code  include  arbitrary  aperture  illumination  and 
incorporating  mode  symmetry.  Computational  run  time  is  reduced  by  nearly  half 
when  symmetry  is  used,  and  significantly  increases  the  applicability  of  this 
method  of  analysis. 

In  order  to  apply  this  analysis  to  multilayer  radomes,  an  equivalent 
surface  impedance  must  be  determined.  The  program  can  be  modified  to  account 
for  radomes  constructed  of  nonuniform  materials  by  assigning  different 
impedances  for  each  segment  of  the  plane  curve  generating  the  BOR. 
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APPENDIX  A.  GAUSS  QUADRATURE 


Computing  the  MM  elements  of  [ll]  and  [z]  (2-19  and  2-28)  as  well  as  the 

gain  requires  the  numerical  evaluation  of  integrals.  MM  requires  numerous 
integration  in  order  to  solve  for  the  current  coefficients  /C  and  l°nJ  of  (2-31).  Not 

only  is  numerical  analysis  necessary  to  perform  the  integrals  in  the  MM.  but  also 
the  matrix  analysis  of  (2-31).  The  large  number  of  separate  integrals  necessitates 
an  integration  technique  which  provides  accuracy  in  a  few  steps  to  avoid 
excessive  run  time.  Gauss  Quadrature  is  a  method  of  numerical  integration  which 
employs  unequally  spaced  intervals  ,  and  provides  accuracy  with  a  few  points. 
[Ref.  4:pp.  154-157] 

Gauss  Quadrature  approximation  of  the  form 

A 

I  =  jf(n)dH,,  (A-l) 


by  the  sum  of  m  terms 


m 


(A-2) 


Jr  =  i 


The  coefficients  Wk  are  weights,  and  the  Hk' s  are  abscissas  determined  from 


”k  = 


b  +  a  b-a c 

- + - 

2  2 


(A-3) 


where  s  are  the  m  number  of  zeros  of  the  m'th  degree  Legendre  Polynomials. 
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APPENDIX  B. 


RADOME.F  LISTING 


The  material  on  the  following  pages  is  a  listing  of  the  program  RADOME.F 
described  in  chapter  III. 
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C*** ***************************************************** **************** 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ROGRAM 

DATE 

LATEST  REVISIOI 
PROGRAMMERS 


r adorns . f  V . 3 

23  January  1992 

31  December  1992 

D.  JENN,  R.  FRAICIS,  K.  KLOPP 


NOTES : 

1.  ARBITRARY  CIRCULAR  APERTURE  ILLUMINATION  SPECIFIED  IN 
FUNCTION  TAPER. 

2.  NEAR  FIELD  COMPONENTS  INCLUDED.  (COMPUTED  IN  SUBROUTINE  GENEX.) 

3.  ANTENNA  IS  AT  THE  BOR  COORDINATE  SYSTEM  ORIGIN. 

4.  LINEARLY  (X-POLARIZED)  ANTENNA. 

5.  CLOSED  FORM  FAR  FIELD  ANTENNA  OPTION  ADDED. 

6.  FLAGS  USED  THROUGHOUT:  (IMP  IS  SET  AUTOMATICALLY) 

IMP=0  perfect  conductor  radome  (lor  test  purposes) 

IPRINT=0  print  pattern  points  to  unit  8 

ISEG=0  print  the  generating  curve  points  to  unit  8 

ICALC=0  compute  current  coefficients  A  pattern 

#0  read  current  coefficients  from  disc  file  given 
by  ’CURCO  and  compute  the  pattern. 

IFF=0  closed  form  far  field  antenna  pattern  is  used  to 
compute  ETF  and  EPF  (SUBROUTINE  ANTFF) 

#0  far  field  computed  using  CIRCRTP  at  distance  ROBS 
ISYM=0  use  mode  symmetry 

#0  turn  off  mode  symmetry  (for  testing) 

7.  THE  FACTOR  ETA  (=377  OHMS)  IS  OMITTED  THROUGHOUT. 

BEGIN  MAIN  PROGRAM********************************************* 


CHARACTER  ITH 

CHARACTER*9  OUT , ETSCATTER , EPSCATTER , CURCO 

CHARACTER* 12  CC 

CHARACTER*6  AN , ETFEED , EPFEED 

CHARACTER*8  GNT , GNPHI , CGP , CPOL , XPOL 

CHARACTER* 14  TPTS ,PPTS .PHIPTS 

COMPLEX  EP,ET,Z( 100000) ,R(1600) ,B(800) ,C(800) ,U 

COMPLEX  UC , ET1 , EP 1 , ET2 , EP2 , EC , EX , ZLO ( 400 ) , ZL ( 2400 ) 

COMPLEX  CIRCR , CIRCT , CIRCP , ETF , EPF 
COMPLEX  EXP1 , EXP2 .CONJG, CEXP , CMPLX , CT(30000) , IMP , JK 
DIMENSION  RH(400) ,ZH(400) ,XT(4) ,AT(4) ,IPS(.800) ,XP(100) ,AP(100) 
DIMENSION  A(100) ,X(100) ,EXP(SOO) ,ANG(500) ,ECP(SOO) 

DIMENSION  EC V ( 500 ) , EXV ( 500 ) ,PHC(500) , PHX(500) 

DIMENSION  ETSCAT(500),EPSCAT(500),ETHF(500),EPHF(500) 

INTEGER  CNPHI, SELECTION 

DATA  IPRINT/O/ , ICALC/1/ , IFF/O/ , ISEG/O/ , ISYM/O/ 

DATA  START, STOP/O. ,360./, PI/3. 1415926/ 

Rad=PI/180. 


ECX=0 . 


BK=2 . *PI 
U=(0. ,1.) 
U0=(0. ,0.) 
UC=-U/4./PI 
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JK=(0. 0,6. 283185307) 

C*****  ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS  *********** 
C*****  (ALL  DATA  FILES  ARE  APPENDED  WITH  THIS  LETTER)  *********** 

C*****  (.m  FILES  ARE  FOR  GENERATING  PLOTS  IN  MATLAB)  *********** 

WRITE(6,*) ’ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS’ 
READ(5 , *)ITH 
OUT= ’ outldbor ’ //ITH 
AN  =’ang’//ITH//’ .m' 

ETFEED= ’etl ’//ITH// 1 .m’ 

EPFEED* * epf ’ //ITH// ’ . m ’ 

ETSCATTER=’etscat’//ITH//’ ,m’ 

EPSCATTER= ’ epscat ’//ITH// ’ .m’ 

CURCO= ’ asciidat ’ //ITH 
CC=’curcoeisdat ’//ITH 
CPOL= ’ etpol ’ //ITH// ’ . m ’ 

XPOL= ’ eppol ’ //ITH// ’ . m ’ 

IF(ICALC.EQ.O)  THEN 

WRITE (6,*) ’ENTER  THE  ANTENNA  RADIUS  (wavelengths) ’ 

READ (5 , *)  ARAD 

C*****SELECT  THE  GEOMETRY  OF  THE  BOR************************************* 
WRITE (6,*) ’SELECT  BOR  GEOMETRY  BY  NUMBER.’ 


WRITE (6,*)’  HUMBER  BOR’ 

WRITE(6,*) ’  1  OGIVE’ 

WRITE(6,*) ’  2  SPHERE’ 

WRITE(6,*) ’  3  CONE’ 

WRITE(6,*) ’  4  DISK’ 

WRITE(6,*) ’  5  PARABOLA’ 

RE AD ( 5 , * ) SELECTI ON 


IF  (SELECTION. EQ.DTHEN 

WRITE(6,*) ’CALLING  SUBROUTINE  FOR  OGIVE’ 

CALL  OGIVE(NP.ZH,RH,BASE,RS,ZP) 

ELSEIF( SELECTION .  EQ .  2 )THEH 

WRITE(6,*) ’CALLING  SUBROUTINE  FOR  SPHERE' 

CALL  TESTSPHERE ( NP , ZH , RH , BASE , RS , ZP ) 

ELSEIF ( SELECTION . EQ . 3 ) THEN 

WRITE(6,*) ’CALLING  SUBROUTINE  FOR  CONE’ 

CALL  CONE ( NP , ZH , RH , B ASE , RS , ZP ) 

ELSEIF ( SELECTION . EQ . 4 ) THEN 

WRITE (6,*) ’CALLING  SUBROUTINE  FOR  DISK’ 

CALL  DISK(NP , ZH , RH , BASE , RS , ZP ) 

ELSLIF (SELECTION .Eq . 5)THEH 

WRITE(6,*) 'CALLING  SUBROUTINE  FOR  PARABOLA’ 

CALL  PARAB ( NP , ZH , RH , DM , FOD , ZP ) 

BASE=DM/2 . 

ENDIF 

IF  (MP.GT.400)THEM 

WRITE(6,*) ’MAXIMUM  NUMBER  OF  POINTS(NP)  IS  400’ 
GOTO  999 
ENDIF 
ENDIF 

C************************************************************ 
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WRITE ( 6 , * ) *  EVTER  TEE  FILEIAMES  gaus###’ 

WRITE(6,*)  ’ . * 

WRITE (6,*) ’EITER  THE  FILEIAME  II  T  (IT)’ 

READ(5,*)GIT 

WRITEC6,*) ’EITER  THE  FILEIAME  II  PHI  (IPHI) ’ 
READ(5,*)GHPHI 

WRITE (6,*) ’ESTER  THE  FILEIAME  II  X  AID  Y’ 

READ(5,*)CGP 

C  OPEH  THE  FILES  FOR  THE  gaus/gaus### 

TPTS=’ gaus/ ’//GIT 

PPTS= 'gaus/ ’//GIPHI 

PHIPTS= ’ gaus/ ’ //CGP 

OPEH ( 1 , FILE=TPTS , STATUS= ’ OLD ’ ) 

OPES ( 2 , FILE=PPTS , STATUS = ’ OLD ’ ) 

OPEH (4 , FILE=PHIPTS , STATUS= ’ OLD ’ ) 

READ(1,*)IT 

IF(HT.GT.4)THEI 

WRITE (6,*) 'MAXIMUM  HUMBER  OF  POIITS(IT)  IS  4’ 

GOTO  999 
EHDIF 

READ (2,*) IPHI 

IF (HPHI . GT . 200 ) THEI 

WRITE (6,*) ’MAXIMUM  HUMBER  OF  POIITS(IPHI)  IS  200’ 
GOTO  999 
EHDIF 

READ(4,*)CHPHI 

IF(CHPHI.GT. 100) THEI 

WRITEC6,*) ’MAXIMUM  HUMBER  OF  POIITS(CIPHI)  IS  100’ 
GOTO  999 
EHDIF 

C  LOAD  THE  WEIGHTS  AID  ABSCISSAS  II  THE  VECTORS. 

DO  1  M=1,IT 

READ( 1 , * , EHD=1 )XT(M) ,AT(M) 

1  COHTIIUE 

DO  2  M=1,HPHI 

READ(2,*,EHD=2)X(M) , A(M) 

2  COHTIIUE 

DO  4  M-l ,CHPHI 

READ(4,*,EHD=4)XP(M) ,AP(M) 

4  COHTIIUE 
CLOSE(l) 

CL0SE(2) 

CLOSE (4) 

WRITE (6,*) ’EITER  THE  PLOTTIIG  IICREMEIT  II  DEGREES’ 
READ(5,*)DT 
IF(ICALC.EQ.O)  THEI 
WRITE (6,*) ’EITER  HIGHEST  MODE’ 

READ(5 , *)  MODES 
HBL0CK=2*M0DES+1 
HR0W=IBL0CK*(2*IP-3 ) 

C  CHECK  FOR  BOUIDS  VIOLATIOI 
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IF(NROW.GT. 30000)  THEM 

WRITE(6.*)  * BOUIDS  VIOLATION  FOR  ARRAY  CT(.)' 

GOTO  999 
END  IF 
ENDIF 

WRITE (6,*) ’ESTER  PHI  (observation)  II  DEGREES' 

READ(5,*)P 

PHI=P*RAD 

IF(ICALC.EQ.O)  THEI 

WRITE ( 6 , * ) ’ EHTER  THE  SCAM  ANGLES  IN  DEGREES:  THETA, PHI' 

READ(5,*)  THD.PHD 

THS=THD*RAD 

PHS=PHD*RAD 

WRITE ( 6, *) 'ENTER  COMPLEX  IMPEDANCE,  OHMS:  (Real , Imag) ’ 

READ ( 5 , * )  IMP 

ELSE 


0*********READ  CURRENT  COEFFICIENTS  IF  ICALC#0************************ 


WRITE(6,*)  'ENTER  LETTER  CODE  FOR  FILE  curcoefsdat’ 
READ(5,*)  ITH 
CC= ’ curcoeisdat ’ //ITH 


OPEN (31 , FILE=CC ) 

WRITE(6,*)  'FILE  OPENED  IS  ’,CC 

READ(31 ,*)  ITH, SELECTION, BASE, NP.MT.MP, I 

RE AD ( 3 1 , * )  NROW , MODES , 1BL0CK , THS , PHS , ARAD , PHI 0 

READ(31 ,*)  IMP , ZP , RS 

READ(31,*)  (RH(L) ,L=1 ,NP) 

READ(31 ,*)  (ZH(L) ,L*l,IP) 

READ (31,*)  (ZLO(L) ,L=1 ,NP) 

READ (31,*)  (CT(L) ,L=1 , NROW) 

CLOSE (31) 

WRITE(6,*) 

ENDIF 


'RUN  CODE  READ  FROM  DISC  :  ’.ITH 


C**********READY  TO  CALCULATE  THE  PATTERN***************************** 
IT=INT ( (STOP-START ) /DT ) + 1 
MHI=M0DES+1 
OPEN ( 8 , FILE=OUT ) 

WRITE ( 8 , 8000 )  2 . *BASE , IP , MODES , RS , ZP 
8000  FORMAT(//, '***  BOR  RADIATION  PATTERN  FOR  A  CIRCULAR’ 

*'  DISC  RADIATOR  USING  GEIEX  ***', 

*//,2X, 'BOR  DIAMETER  (WAVELENGTHS)*’ ,F5. 2,/, 2X, 

*  'NUMBER  OF  GENERATING  POINTS  (NP)=’ ,I4,/,2X, 

*  'HIGHEST  MODE  NUMBER  USED  (MODES)*’ ,13, 

* // , 2X , ’ SURFACE  RADIUS’ ,F5.2,/ ,2X, ’ ZPRIME ’ , FS . 2) 

WRITE (8, 30)  NT.NPHI .CNPHI 
30  F0RMAT(/,12X, ’  NT  NPHI  CNPHI’, 

*  /,10X,I5,2X,I5,2X,I5) 

IF(ISEG.EQ.O)  WRITE(8 , 1300) 

1300  FORMAT(/, 10X, 'INDEX ’ ,8X, ’Z(I) ’ , 10X , ’RHO(I) ' , 12X, 'SURF  IMPED’) 
MP=NP-1 
MT=MP-1 
N=MT+MP 
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DO  52  1=1, IP 

IF(ABS(ZH(I)).LT. .001)  ZH(I)=0. 

IF(ABS(RH(I)).LT. .001)  RH(I)=0. 

ZHB=ZH(I)/BK 

RHB=RH(I)/BK 

C  ASSIGH  SURFACE  IMPEDAICE  AT  THIS  POUT.  THE  SURFACE  IMPED AICE  OF  SEGMEIT 
C  I  IS  ZLO(I) 

IF( ICALC . Eq . 0)  Z10( I ) =IMP/ ( 120 . *PI ) 

IF(ISEG.EQ.O)  WRITE(8 ,8004)  I.ZHB.RHB.IMP 
52  COITIIUE 

8004  FORMAT ( 1 IX , 14 , 4X , F8 . 3 , 8X , F8 . 3 , 6X , 2F8 . 2 ) 

C************************************************************************ 

C  BIG  LOOP  *>***************************** 

C 

C  MODE  LOOP  TO  CALCULATE  THE  CURREIT  COEFFICIESTS.  POS  AID  IEG  MODES  * 

-  :OHE  II  THE  SAME  ITERATIOI  OF  THE  LOOP 


C  T 

C  400  COHTIHUE  <  **************************** 

C**********************************************ZLO AD, ZMAT,GEIEX, DECOMP, SOLVE 
IF ( ICALC. Eq.O)  THEI 

IF(CABS(IMP) .IE.O)  CALL  ZLOAD(IP,RH,ZH,ZLO,ZL) 

DO  400  M=1,MHI 
HM=M-1 

CALL  ZMAT ( HM , IN , IP , IPHI , IT , RH , ZH , X , A , XT , AT , Z ) 


C*******************************************************POSITIVE  MODE  +n 


C 

C 

C 

c 

c 

c 

c 


ZMAT  = 


Z 

+n 

phi.t 

Z 

+n 


t  ,phi 
Z 

+n 

phi, phi 
Z 

+n 


C*********************************************************************** 


c  MODE  SYMMETRY  II  n  FOR  ZMAT  EXISTS  BUT  IS  IOT  USED  II  THIS  VERSIOI. 

C  MODE  SYMMETRY  II  I  USED  DIRECTLY  FOR  PHIS  OF  0,90,180,  AID  270. 

C*********************************************************************** 


C 

1 

t,t 

t.phi  1 

1  t,t 

t.phi 

c 

1 

Z 

Z  1 

1  z 

-Z 

c 

1 

-n 

-a  1 

1  +n 

+n 

c 

1 

|  = 

1 

c 

1 

phi.t 

phi, phi  ! 

1  phi , t 

phi, phi 

c 

1 

Z 

Z  1 

1  -Z 

Z 

c 

1 

-n 

-n  1 

1  +n 

+n 

C* ******************************** *******************CURREIT  COEFFICIEITS 

C  _ -1 


C 

1 

1  1 

1 

1  1 

C 

! 

I  I 

1  Z  1 

0  1 

0 

1  1 

V 

C 

1 

+nl 

1  +n  1 

1 

+R 

C 

1 

1 

1  1 

1 

_ 

_l 

C 

1 

1  1 

1 

1  I 

C 

CT  =  1 

I  1  = 

1  0  I 

Z  1 

0 

1  * 

V 
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0 


0 


c 

c 

c 

c 

c 

c 


I 


I 


0 


0 


0 


z 


V 


-n 


C***************************************************POSITIVE  MODE  +n 

IF(CABS(IMP) .IE.O)  CALL  ZTOT(MT,KP,ZL,Z) 

CALL  GEIEX(IM , HP , IT , IPHI , CIPHI , XT , AT , X , A , 

*  XP,AP,THS,PHS,ARAD,RH,ZH,B) 

CALL  DECOMP(I,IPS,Z) 

CALL  SOLVE (H , IPS , Z . B , C ) 

ITOPl=MODES-IM 
jrT0P2=MBL0CK-  (  ITOP 1 + 1 ) 

IS2=IT0P1*I 

IS1=IT0P2*I 

NRQW=IBLOCK*I 


C******************************** ************************ ******* 

C  STORE  CURREIT  COEFFICIESTS  II  OIE  LOIG  VECTOR  CT( . ) .  MOST  IEGATIVE 
C  MODE  IS  AT  THE  TOP;  MOST  POSITIVE  AT  THE  BOTTOM.  UIIT  30  IS  FORMATED 
C  OUTPUT;  UIIT  31  IS  FREE  FORMATED  (USED  BY  GAII  PROGRAM). 
C*************************************************************** 

OPEI ( 30 , FILE=CURCO ) 

WRITE (30, 5027) ’Mode( 1 ,IM, ’ )  =  '.IM.’IM  =  \IM 
DO  401  L-l ,1 

WRITE (30, 5026)  >CT( ’ ,L+  IS1 , ’)  *  *,C(L) 

401  CT(L+IS1)=C(L) 

5027  FORMAT (/A, I, A, I/A, I/) 

5026  FORMAT (5X ,  A ,  14 ,  A ,  F10 . 4 ,  ’  +  ( ’ ,F10.4, ’ )*i’ ) 

C* ************* *EID  OF  POSITIVE  MODE*******BEGII  IEGATIVE  MODE  -n 
C  MODE  SYMMETRY  IS  USED  IF  ISYM=0 
IF(IM.IE.O)  THEI 
IMI=-IM 

C  USE  BRUTE  FORCE  IF  ISYM  #  0 
IF(ISYM.IE.O)  THEI 

CALL  GEIEX(IMI , IP ,IT , IPHI , CIPHI , XT , AT , X , A , 

*  XP , AP , THS , PHS , ARAD , RH , ZH , B ) 

CALL  ZMAT ( IMI , IMI , IP , IPHI , IT , RH , ZH , X , A , XT , AT , Z ) 

IF(CABS(IMP).IE.O)  CALL  ZTOT(MT,MP,ZL,Z) 

CALL  DECOMP (I, IPS, Z) 

CALL  SOLVE(I,IPS,Z,B,C) 

EIDIF 

WRITE ( 30 , 5027 ) 1  Mode ( ’ , IMI , ’ )  =  »,IMI,’IM  =  >,IM 

DO  402  L=1 ,MT 

CT(L+IS2)=+C(L) 

WRITE (30, 5026)  *CT( ' .L+IS2 , ’ )  =  \CT(L+IS2) 

402  COITIIUE 

DO  403  L=MT+1,I 
IF(ISYM.IE.O)  CT(L+IS2)=+C(L) 

IF(ISYM.EQ.O)  CT(L+IS2)=-C(L) 

WRITE (30, 5026)  *CT( ’ ,L+IS2 , ’ )  =  ’ ,CT(L+IS2) 
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403  CONTINUE 

0*****************************************END  OF  NEGATIVE  MODE  -n 
EBDIF 

400  CONTINUE 
CL0SE(30) 

C*<  :i******************************************************************** 

C  WRITE  THE  VECTOR  OF  CURREHT  COEFFICIEHTS  01  DISC  FOR  GAII  CALC. 

C  (FREE  FORMAT  TO  UNIT  31.) 

C*********************************************************************** 

0PEM(31,FIL£=CC) 

WRITE (31,*)  ITH , SELECTION , B ASE , IP , MT , MP , B 
WRITE (31,*)  BROW , MODES , NBLOCK , THS , PHS , ARAD , PHI 
WRITE(31,*)  IMP , ZP , RS , IFF 
WRITE(31,*)  (RH(L) ,L=1 ,BP) 

WRITE(31 , *)  (ZH(L) ,L=1 ,BP) 

WRITE(31,*)  (ZLO(L) ,L=1 , BP) 

WRITE(31,*)  (CT(L) ,L=1 .BROW) 

WRITE(31 ,*)  GBT.GBPHI.CGP 

WRITE(6,*)  'CURRENT  COEFFICIEHTS  WRITTEN  ON  DISC’ 

CL0SE(31) 

END  IF 

C**********************************************************«************* 
C  END  OF  BIG  LOOP  *>************************* 
C  *  v 

C*******<  ************  *************** 

c*** *************************************************** ****************** 
C 

C  BEGIN  PATTERN  LOOP  FROM  START  TO  STOP  IN  INCREMENTS  OF  DT  (ALL  IB  DEG) 

C 

DO  500  1=1, IT 

THETA=START  +  DT*FL0AT(I-1) 

THX=THETA*RAD 

PHR=PHI 

IF (THETA . GT . 180 . )  THEN 
PHR=PHI+PI 

THX= ( 360 . -THET A ) *RAD 
ENDIF 

ET1=(0 . ,0.) 

EP1=(0. ,0. ) 

ET2=(0 . ,0. ) 

EP2=(0. ,0.) 

DO  300  H=1,MHI 
NM=M-1 

EXP 1 =CEXP ( CMPLX ( 0 . ,FLOAT(NM)*PHR) ) 

EXP2=C0B  J  G ( EXP 1 ) 

C*** ************************ ********************************PLABE 

CALL  PLANE(BM,NM,NP,BT,RH,ZH,XT, AT,TBX,R) 

NTOP 1 =M0DES-NM 

NT0P2 = NBLOCK- ( NTOP 1 + 1 ) 

BS2=NT0P1*B 

BSi=BT0P2*B 
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DO  250  L=1,MT 

ET1=ET1*R(L)*CT(L+IS1)*EXP1 
EP 1 =EP 1 +R ( L+I ) *CT( L+IS 1 ) *EXP 1 
IF(IM.EQ.O)  GOTO  250 
ET1=ET1+R(L)*CT(L+IS2)*EXP2 
EP1=EP1-R(L+I)*CT(L+IS2)*EXP2 
250  COITIHUE 

DO  260  L=1 ,MP 

ET2=ET2+R(L+MT)*CT(L+IS1+MT)*EXP1 
EP2=EP2+R(L+MT+M): CT(L+IS1+MT) *EXP1 
IF(HM.EQ.O)  GO  TO  260 
ET2=ET2-R(L+MT)*CT(L+HS2+MT)*EXP2 
EP2=EP2+R(L+MT+I)*CT(L+NS2+MT)*EXP2 
260  COITIHUE 
300  COITIHUE 

C  FEED  COHTRIBUTIOH  II  THE  FAR  FIELD  IS  ETF.EPF 

C*** ********** ***************************************BESSJ1 . 

C  USE  CLOSED  FORM  FEED  EXPRESSIOI  FROM  SUBROUTIIE  AITFF  IF 
C  IFF=0 ;  OTHERWISE  USE  BRUTE  FORCE  EVALUATIOH  FROM  CIRCRTP 
C 

ROBS- 1000. 

IF(IFF.EQ.O)  THEI 

CALL  AITFF (THX.PHR, THS , PHS , ARAD .ETF.EPF) 

ELSE 

RHB=ROBS*SII (THX ) 

ZHB=ROBS*COS (THX ) 

CALL  CIRCRTP (CIPHI , XP , AP . ARAD , THS . PHS . 

*  PHR , RHB , ZHB , CIRCR , CIRCT , CIRCP ) 

C  REMOVE  THE  1/R  DEPEIDEICE  BECAUSE  EXP(-jkR)/R  IS  OMITTED  II 
C  THE  SCATTERED  FIELDS  ET  AID  EP 
ETF=CIRCT*ROBS 
EPF=CIRCP*ROBS 
EHDIF 

C**************************************************** 

ETHF (I)=CABS( ETF ) 

EPHF ( I ) =CABS (EPF ) 

ET=(ET1+ET2)*UC 
EP=(EP1+EP2)*UC 
ETSCAT ( I ) =C ABS (ET) 

EPSCAT ( I ) =C ABS ( EP) 

C  TOTAL  E-THETA  AID  E-PHI  COMPOIEITS 
EC=ET+ETF 
EX=EP+EPF 
ECV ( I ) =CABS (EC ) 

EXV(I)=CABS(EX) 

ECR-REAL(EC) 

ECI=AIMAG(EC) 

EXR=REAL(EX) 

EXI=AIMAG(EX) 

PHC(I)=ATAI2(ECI ,ECR+1 . «-10)/RAD 
PHX(I)=ATAI2(EXI,EXR+1.«-10)/RAD 
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ANG(I)=THETA 

ECX= AMAX1 ( ECX , ECV ( I ) ,EXV(I)) 

500  CONTINUE 

WRITE (6,*)  ’MAX  E  VALUE=',ECX 
WRITE (8,103)  P , THS/RAD . PHS/RAD , ARAD , ECX 
103  FORMAT (/, 10X, ’PHI  OF  RECEIVER  (DEG)= ’ ,F8. 2,/ , 10X, 

*  ’ ANTENNA  SCAM:  THETA  (DEG)= ’ ,F8 . 2 ,/ , 10X, 

*  ’  PHI  (DEG)=* ,F8.2,/,10X, 

*  ’ ANTENNA  RADIUS:  ARAD=* ,F8.3,/,10X, 

*  ’MAXIMUM  FIELD  VALUE  (V/M)= ’ .E15.5) 

IF(IFF.EQ.O)  THEM 

WRITE(8 , 113) 

ELSE 

WRITE(8,213)  ROBS 
EHDIF 

113  F0RMAT(/,10X, ’CLOSED  FORM  FAR  FIELD  PATTERN  FROM  AITFF  IS  USED’) 
213  FORMAT ( / , 10X , ’ FAR  FIELD  COMPUTED  USIHG  CIRCRTP,  ROBS=’ ,F10.2) 

C  STORE  DATA  FOR  NORMALIZED  CO-  AID  CROS-POLARIZED  PATTERNS 
DO  600  1=1, IT 
ECP(I)=(ECV(I)/ECX)**2 
EXP  ( I )  =  (EXV  ( I )  /ECX )  **2 
ECP(I)=AMAX1(ECP(I) .l.E-6) 

EXP(I)=AMAX1(EXP(I) , 1 .E-6) 

ECP(I)=10.*AL0G10(ECP(I) ) 

EXP(I)=10.*ALOG10(EXP(I) ) 

600  CONTINUE 

IF(IPRINT.EQ.O)  THEN 
WRITE(8 ,5015) 

5015  FORMAT ( // , 7X , ’ ANGLE ’, 15X ,’ CO-POL ’ , 25X , ’ X-POL ’ , / , 7X , 

1 ’ (DEG) ’ ,4X,2( ’ (VOLTS) ’ ,4X, ’ (DEG) ’ ,3X, ’ (DB-REL) ’ ,4X)) 

DO  9000  L=1 , IT 

WRITE(8 ,5016)  ANG(L) ,ECV(L) ,PHC(L) ,ECP(L) ,EXV(L) ,PHX(L) 

1 , EXP(L) 

5016  FORMAT (5X,F6.2,3X,2(F8.4,3X,F7.2,3X,F7.2,3X)) 

9000  CONTINUE 

ENDIF 

0PEN(2 ,lil«=AN) 

OPEN ( 3 , 1 ile=CPOL ) 

0PEN(4,f ile=XPOL) 

OPEN (7 , fiI«=ETSCATTER) 

OPEN (8 , 1 ile=EPSCATTER) 

0PEN(9,file=ETFEED) 

OPEN ( 1 0 , f il e=EPFEED ) 

DO  9097  1=1, IT 
WRITE(2, 5019)1,  ANG(I) 

WRITE (3, 5020)1,  ECP(I) 

WRITE(4, 5021)1,  EXP(I) 

WRITE (7,5022)1,  ETSCAT(I) 

WRITE(8, 5023)1,  EPSCAT(I) 

WRITE(9 ,5024)1 ,  ETHF(I) 

9097  WRITE( 10,6025) I.  EPHF(I) 


5019  FORMAT ( ’  ANG( ’ , I , ’ ,F8. 3, ’ ; ’) 

5020  FORMAT ( ’ECP(* ,1, ’)=’ ,F8.3, ’ ;  ■) 

5021  FORMAT ( ’EXPO  ,1 ,  '  ,F8. 3 ,  ’ ;  ’ ) 

5022  FORMAT ( ’ ETSCAT (’,I,’)=’,F8.3, * ; ') 

5023  FORMAT (’EPSC AT (’ ,1, ’)=’,F8.3, *) 

5024  FORMAT(’ETHF(\I,  ’)  =  \F8.3,  *; ') 

5025  FORMAT (’EPHFC ,1, ’)=’ ,F8.3, ’) 

CL0SE(2) 

CL0SE(3) 

CL0SE(4) 

CL0SE(7) 

CL0SE(8) 

CLOSE (9) 

CLOSE (10) 

999  CONTINUE 
STOP 
F.HD 

0********************=**********************************END  OF  MAII  PROGRAM. 

C 

C********************************************************SUBROUTIFE  ZMAT. 

C  REFERENCE :  AN  IMPROVED  E-FIELD  SOLUTION  FOR  CONDUCTING  BOR 
C  J.R.MAUTZ  AND  R.F.  HARRINGTON 

C  TECHNICAL  REPORT  TR-80-1 

C  ROME  AIR  DEVELOPMENT  CENTER 

C  CONTRACT  NO.  F-30602-79-C-0011 

SUBROUTINE  ZMAT (Ml . M2 , NP , NPHI . NT , RH , ZH , X , A , XT , AT , Z ) 

C 

C  COMPUTE  THE  MM  IMPEDANCE  MATRIX  ELEMENTS.  THIS  IS  FROM  MAUTZ  AND 
C  HARRINGTON  (NO  CHANGES). 

C 

COMPLEX  Z( 100000) ,U1 ,U2,U3,U4,U5,U6,U7,U8,U9,UA,UB,G4A(4) ,G5A(4) 
COMPLEX  G6 A ( 4 ) , G4B ( 4 ) , G5B ( 4 ) , G6B (4 ) , H4A , H5 A , H6A , H4B , H5B 
COMPLEX  H6B,UC,UD,GA(100) ,GB(100) 

DIMENSION  RH(400) ,ZH(400) ,X(100) ,A(100) ,XT(4) ,AT(4) 

DIMENSION  RS(400) ,ZS(400) ,D(400) ,DR(400) ,DZ(400) 

DIMENSION  DM(400) ,C2(100) ,C3(100) ,R2(4) ,Z2(4) 

DIMENSION  C4(100) ,C5(100) ,C6(100) ,Z7(4) ,R7(4) ,Z8(4) ,R8(4) 

CT=2. 

CP=.l 

DO  10  1=2 ,NP 
12=1-1 

RS(I2)=.5*(RH(I)+RH(I2)) 

ZS(I2)= . 5*(ZH(I)+ZH(I2) ) 

Dl= . 5*(RH(I)-RH(I2) ) 

D2=. 5*(ZH(I)-ZH(I2) ) 

D(I2)=SQRT(D1*D1+D2*D2) 

DR(I2)=D1 

DZ(I2)=D2 

IF(RS(I2) . EQ.O.  )  RS(I2)  =  1 . 

DM(I2)=D(I2)/RS(I2) 

10  CONTINUE 


75 


M3=M2-M1+1 

H4=M1-1 

PI2=1. 570796 

DO  11  K=1 ,HPHI 

PH=PI2*(X(K)+1 . ) 

C2(K)=PH*PH 

SH=SIH(.5*PH) 

C3(K)=4.*SN*SK 

A1=PI2*A(K) 

D4=.5*A1*C3(K) 

D5=A1*C0S(PH) 

D6=A1*SII(PH) 

M5=K 

DO  29  M=1 ,M3 
PHM=(M4+M)*PH 
A2=C0S(PHM) 
C4(M5)=D4*A2 
C5(M5)=D5*A2 
C6(M5)=D6*SII(PHM) 
M5=M5+IPHI 
29  COITIHUE 
11  COHTIIUE 
MP=IP  1 
MT=HP-1 
H=MT+MP 
I2I=KT*I 
I2=f*« 

Ul=(0.,.5) 

U2=(0. ,2.) 

JI=-1-I 

DO  15  jq=l,MP 

Kq=2 

iF(jq.Eq.i)  Kq=i 
iF(jq.Eq.MP)  Kq=3 
Rl=RS(Jq) 

zi=zs(jq) 

Di=D(jq) 

D2=DR( jq) 

D3=DZ(jq) 

D4=D2/R1 
D5=DM( jq) 

SV=D2/D1 

CV=D3/D1 

T6=CT*D1 

T62=T6*D1 

T62=T62*T62 

R6=CP*R1 

R62=R6*R6 

DO  12  L=1 ,IT 

R2(L)=R1+D2*XT(L) 

Z2(L)=Z1+D3*XT(L) 


12  COHTIHUE 
U3=D2*U1 
U4=D3*U1 
DO  16  IP=1 ,MP 
R3=RS(IP) 

Z3=ZS(IP) 

R4=R1-R3 

Z4=Z1-Z3 

FM=R4*SV+Z4*CV 

PHM=ABS(FM) 

PH=ABS(R4*CV-Z4*SV) 

D6=PH 

IF(PHM.LE.Dl)  GO  TO  26 
D6=PHM-D1 

D6=SQRT(D6*D6+PH*PH) 

26  IF(IP.Eq.JQ)  GO  TO  27 
KP=1 

IF(T6.GT.D6)  KP=2 
IFCR6.GT.D6)  KP=3 
GO  TO  28 

27  KP=4 

28  GO  TO  (41,42,41,42) ,KP 
42  DO  40  L=1 ,IT 

D7=R2(L)-R3 

D8=Z2(L)-Z3 

Z7(L)=D7*D7+D8*D8 

R7(L)=R3*R2(L) 

Z8(L)=.25*Z7(L) 

R8(L)=.25*R7(L) 

40  COHTIHUE 

Z4=R4*R4+Z4*Z4 
R4=R3*R1 
R5=.5*R3*SV 
DO  33  K=1 ,IPHI 
A1=C3(K) 

RR=Z4+R4*A1 
UA=0 . 

UB=0 . 

IFCRR.LT.T82)  GO  TO  34 
DO  35  L=1 ,IT 
R=SQRT(Z7  (L)'*-R7(L)*A1) 
SI=-SI1T(R) 

CS=COS(R) 

UC=AT(L)/R*CH?LX(CS,S*) 
UA=UA+UC 
UB=XT(L)*UC+UB 
35  COITIIUE 
GO  TO  36 
34  DO  37  L=1 ,IT 

R=SqRT(Z8(L)-t-R8(L)*Al) 

S«=-SI1(R) 


CS=COS(R) 

UC=AT ( L ) /R*SR*CMPLX ( -SR , CS ) 

UA=UA+UC 

UB=XT(L)*UC+UB 

37  COHTIRUE 
A2=FM+R5*A1 
D9=RR-A2*A2 
R=ABS(A2) 

D7=R-D1 

D8=R+D1 

D6=SQRT(D8*D8+D9) 

R=SQRT(D7*D7+D9) 

IF(D7.GE.O. )  GO  TO  38 
A l=ALOG ( (D8+D6 ) * (-D7+R) /D9 ) /D 1 
GO  TO  39 

38  Al=ALOG( (D8+D6)/(D7+R) )/Dl 

39  UA=A1+UA 

UB=A2*(4./(D6+R)-A1)/D1+UB 
36  GA(K)=UA 
GB(K)=UB 
33  CORTIHUE 
K1=0 

DO  45  M=1 ,M3 
H4A=0 . 

H5A=0. 

H6A=0. 

H4B=0 . 

H5B=0 . 

H6B=0 . 

DO  46  K=1,HPHI 

K1=K1+1 

D6=C4(K1) 

D7=C5(K1) 

D8=C6(K1) 

UA=GA(K) 

UB=GB(K) 

H4A=D6*UA+H4A 
H5A=D7*UA+H5A 
H6A=D8*UA+H6A 
H4B=D6*UB+H4B 
HSB=D7*UB+H5B 
H6B=D8*UB+H6B 
46  corriiUE 
G4A(M)=H4A 
G5A(M)=H5A 
G6A(M)=H6A 
G4B(H)=H4B 
G5B(M)=H5B 
G6B(M)=H6B 
45  CORTIIUE 

IF(KP .RE.4)  GO  TO  47 
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A2=D1/(PI2*R1) 

D6=2./D1 
D8=0 . 

DO  63  K=1 ,NPHI 
A1=R4*C2(K) 

R=R4*C3(K) 

IFCR.LT.T62)  GO  TO  64 
D7=0  . 

DO  65  L=1 ,NT 

D7=D7+AT(L)/SQRT(Z7(L)+A1) 

65  CONTINUE 
GO  TO  66 

64  A1=A2/(X(K)+1 . ) 

D7=D6*AL0G(A1+SQRT(1 . +A1*A1) ) 

66  D8=D8+A(K)*D7 
63  CONTINUE 

A1=.5*A2 

A2=1./A1 

D8=-PI2*D8+2./Rl*(BLOG(A2)+A2*BLOG(Al)) 

DO  67  M=1,M3 
G5A(M)=D8+G5A(H) 

67  CONTINUE 
GO  TO  47 

41  DO  25  M-l ,M3 
G4A(M)=0. 

G5A(M)=0. 

G6A(M)=0. 

G4B(M)=0. 

G5B(M)=0. 

G6B(M)=0. 

25  CONTINUE 
DO  13  L=1 ,NT 
A1=R2(L) 

R4=A1-R3 

Z4=Z2(L)-Z3 

Z4=R4*R4+Z4*Z4 

R4=R3*A1 

DO  17  K=1,NPHI 

R=SQRT(Z4+R4*C3(K) ) 

SN=-SIH(R) 

CS=COS(R) 

GA(K)=CMPLX(CS,SN)/R 
17  CONTINUE 
D6=0 . 

IF(R62.LE.Z4)  GO  TO  51 
DO  62  K=1 ,NPHI 
D6=D6+A(K)/SQRT(Z4+R4*C2(K) ) 

62  CONTINUE 

Z4=3. 141593/SQRT(Z4/R4) 

D6=-PI2*D6+AL0G(Z4+SQRT( 1 . +Z4*Z4) )/SQRT(R4) 
51  A1=AT(L) 
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A2=XT(L)*A1 

K1=0 

DO  30  M=1 ,M3 
U5=0 . 

U6=0. 

U7=0. 

DO  32  K=1 ,NPHI 
UA=GA(K) 

K1=K1+1 

U5=C4(K1)*UA+U5 
U6=CS(K1)*UA+U6 
U7 =C6 (K 1 ) *UA+U7 
32  CONTINUE 
U6=D6+U6 

G4A(M)=A1*U5+G4A(M) 
G5A(M)=A1*U6+G5A(M) 
G6A(M)=A1*U7+G6A(M) 
G4B (M ) =A2*U5+G4B (M ) 
G5B(M)=A2*U6+G5B(M) 
G6B(M)=A2*U7+G6B(M) 
30  CONTINUE 
13  CONTINUE 
47  A1=DR(IP) 

UA=A1*U3 

UB=DZ(IP)*U4 

A2=D(IP) 

D6=-A2*D2 

D7=D1+A1 

D8=D1*A2 

JM=JN 

DO  31  M=1 ,M3 

FM=M4+M 

A1=FM*DM(IP) 

H5A=G5A(M) 

H5B=G5B(M) 

H4A=G4A(M)+H5A 

H4B=G4B(N)+H5B 

H6A=G6A(M) 

H6B=G6B(M) 

U7=UA*HSA+UB*H4A 

U8=UA*H5B+UB*H4B 

U5=U7-U8 

U6=U7+U8 

U7=-U1*H4A 

U8=D6*H6A 

U9=D6*H6B-A1*H4A 

UC=D7*(H6A+D4*H6B) 

UD=FM*D6*H4A 

K1=IP+JM 

K2*K1+1 

K3=K1+N 
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K4=K2+I 

KS=K2+MT 

K6=K4+MT 

K7=K3+H2H 

K8=K4+H2H 

GO  TO  (18, 20, 19), KQ 

18  Z(K6)=U8+U9 
IF(IP.EQ.l)  GO  TO  21 
Z (K3 ) =Z (K3 ) +U6-U7 

Z (K7 ) =Z (K7 ) +UC-UD 
IF(IP.EQ.MP)  GO  TO  22 

21  Z(K4)=U6+U7 
Z(K8)=UC+UD 
GO  TO  22 

19  Z(K5)=Z(K5)+U8-U9 
IF(IP.Eq.l)  GO  TG  23 
Z(Kl)=Z(Kl)+U5+U7 
Z(K7)=Z(K7) +UC-UD 
IF(IP.EQ.MP)  GO  TO  22 

23  Z (K2 ) =Z (K2 ) +U5-U7 
Z(K8)=UC+UD 

GO  TO  22 

20  Z(K5)=Z(K5)+U8-U9 
Z(K6)=U8+U9 
IF(IP.Eq.l)  GO  TO  24 
Z(K1)=Z(K1)+U5+U7 
Z(K3)=Z(K3) +U6-U7 

Z (K7 ) =Z (K7 ) +UC-UD 
IF(IP.Eq.MP)  GO  TO  22 

24  Z(K2 ) =Z (K2 ) +U5-U7 
Z(K4)=U6+U7 
Z(K8)=UC+UD 

22  ZCK8+MT) =U2* (D8* (H5A+D4*H5B ) -A1*UD) 

JM=JM+H2 

31  COSTISUE 
16  COITIHUE 
JI=JI+I 
15  COITIIUE 
RETURI 
ElfD 

C********************************************************SUBR0UTHE  SOLVE. 
C  REFEREICE  :  TECHIICAL  REPORT  TR-80-1 
SUBROUTIIE  SOLVE(l, IPS ,UL,B ,X) 

C 

C  SEE  MAUTZ  *  HARRIIGTOI  FOR  DETAILS 
C 

COMPLEX  UL(IOOOOO) ,B(800) ,X(800) ,SUM 
DIMEHSIOI  IPS(800) 

*P=I+1 

IP=IPS(1) 

X(1)=B(IP) 
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DO  2  1=2,1 
IP=IPS(I) 

IPB=IP 
IM1=I-1 
SUM=(0. ,0.) 

DO  1  J=1 , IM1 
SUH=SUM+UL( IP) *X( J ) 

1  IP=IP+H 

2  X(I)=B(IPB)-SUM 
K2=I*(H-1) 

IP=IPS(I)+K2 

X(I)=X(I)/UL(IP) 

DO  4  IBACK=2,N 

I=HP-IBACK 

K2=K2-I 

IPI=IPS(I)+K2 

IP1=I+1 

SUM=(0 . .0.) 

IP=IPI 

DO  3  J=IP1,H 
IP=IP+I 

3  SUM=SUM+UL(IP)*X(J) 

4  X(I)=(X(I)-SUM)/UL(IPI) 

RETURI 

EHD 

C********************************************^*SUBROUTIIE  DECOMP 
C  REFEREICE:  TECHIICAL  REPORT  TR-80-1 
SUBROUTIIE  DECOMP(I,IPS,UL) 

C 

C  SEE  MAUTZ  k  HARRIHGTON  FOR  DETAILS 
C 

COMPLEX  UL(IOOOOO), PIVOT, EM 
DIMEISIOI  SCL(800),IPS(800) 

DO  5  1=1,1 
IPS(I)=I 
RI=0 . 

J1=I 

DO  2  J=1,I 

ULM= ABS (REAL(UL( J1 ) ) ) +ABS ( AIMAG (UL( J 1 ) ) ) 

J1=J1+* 

IF(RI-ULM)  1,2,2 

1  RI=ULM 

2  COITIIUE 
SCL(I)=1./RI 

5  COITIIUE 
IM1=I-1 
K2=0 

DO  17  K=1,IM1 
BIG=0 . 

DO  11  I=K,I 
IP=IPS(I) 
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IPK=IP+K2 

SIZE=(ABS(REAL(UL(IPK)))+ABS(AIMAG(UL(IPK))))*SCL(IP) 

IF(SIZE-BIG)  11,11,10 

10  BIG=SIZE 
IPV=I 

11  COHTIHUE 
IF(IPV-K)  14,15,14 

14  J=IPS(K) 

IPS(K)=IPS(IPV) 

IPS(IPV)=J 

15  KPP=IPS(K)+K2 
PIVOT=UL(KPP) 

KP1=K+1 

DO  16  I=KP1,H 
KP=KPP 
IP=IPS(I)+K2 
EM=-UL( IP) /PIVOT 
18  UL(IP)=-EM 
DO  16  J=KP1,H 
IP=IP+M 
KP=KP+* 

UL(IP)=UL(IP)+EM*UL(KP) 

16  COHTIHUE 
K2=K2+M 

17  COHTIHUE 
RETURH 
EHD 

O**********************************«*************FUHCTI0H  BLOG 
C  REFEREHCE:  TECHHICAL  REPORT  TR-80-1  ; (page  56) 

FUHCTIOH  BLOG(X) 

IFtt.GT. .1)  GO  TO  1 
X2=X*X 

BLOG= ( ( . 075*X2- . 1666667 ) *X2+ 1 . ) *X 
RETURH 

1  BL0G=AL0G(X+SQRT(1 ,+X*X) ) 

RETURH 

EHD 

C* *********************************************** ♦♦♦♦SUBROUTIHE  PLAHE 
C  REFEREHCE:  TECHHICAL  REPORT  TR-80-1  ; (pages  57-64) 

SUBROUTIHE  PLAHE ( Ml , M2 , HP , HT , RH , ZH , XT , AT , THR , R) 

C 

C  PLAHE  HAVE  EXCITATIOH  VECTOR  IH  THE  FAR  FIELD.  FROM  MAUTZ  AID  HARRIIGTOI. 
C  MODIFIED  TO  DO  OILY  OIE  AIGLE  AID  FREQUEICY  PER  CALL. 

C  **  MEW  BESSEL  FUHCTIOH  ROUTIIE  FROM  RUM.  RECIPES  ** 

C 

COMPLEX  R(1600) , U,U1 ,UA,UB ,FA( 1500) ,FB(1500) ,F2A,F2B,F1A,F1B 
COMPLEX  U2.U3.U4, US 

DIMEHSIGH  RH(400) ,ZH(400) ,XT(4) ,AT(4) ,R2(4) ,Z2(4) 

MP=HP-1 

MT=MP-1 

H*MT+MP 
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12=2*1 

CC=COS(THR) 

SS=SII(THR) 

U=(0. ,1.) 

Ul=3. 141593*U**M1 

M3=M1+1 

M4=M2+3 

IF(Ml.EQ.O)  M3=2 

M5=Ml+2 

M6=M2+2 

DO  12  IP=1,MP 

K2=IP 

I=IP+1 

DR=.5*(RH(I)-RH(IP) ) 
DZ=.5*(ZH(I)-ZH(IP) ) 
D1=SQRT(DR*DR+DZ*DZ) 

Rl=. 25*(RH(I)+RH(IP) ) 

IFCRl.Eq.O.)  Rl=l . 
Z1=.5*(ZH(I)+ZH(IP)) 

DR=.5*DR 

D2=DR/R1 

DO  13  L=1 ,IT 

R2(L)=R1+DR*XT(L) 

Z2(L)=Z1+DZ*XT(L) 

13  COITIIUE 
D3=DR*CC 
D4=-DZ*SS 
D5=D1*CC 
DO  23  M=M3,M4 
FA(M)=0. 

FB(M)=0. 

23  COITIIUE 
DO  15  L=1 ,IT 
X=SS*R2(L)*2 . 

ARG=Z2(L)*CC 

UA=AT(L)*CMPLX(COS(ARG) , SII(ARG) ) 
UB=XT(L)*UA 

C  THIS  LIRE  REPLACES  HARRIIGTOI’S  BLOCK 
DO  25  M=H3,M4 
BES=BESS J (M-2 , X ) 
FA(M)=BES*UA+FA(M) 
FB(M)=BES*UB+FB(M) 

25  COITIIUE 
15  COITIIUE 

IF(Ml.IE.O)  GO  TO  26 
FA(1)=-FA(3) 

FB(1)=-FB(3) 

26  UA=U1 

DO  27  N=K5,N6 

M7=M-1 

H8=M+1 
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F2A=UA*(FA(M8)+FA(M7) ) 

F2B=UA*(FB(M8)+FB(M7) ) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7) ) 

F1B=UB*(FB(M8)-FB(M7) ) 

U4*D4*UA 

U2=D3*F1A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+H 

K5=K2+I 

R(K2+MT ) =-D5* (F2A+D2*F2B ) 

R(K5+MT)=D1*(F1A+D2*F1B) 

IF(IP.EQ.l)  GO  TO  21 
R(K1)=R(K1)+U2-U3 
R(K4 ) =R (K4 ) +U4-U5 
IF(IP.EQ.MP)  GO  TO  22 

21  R(K2)=U2+U3 
R(K5)=U4+U5 

22  K2=K2+H2 
UA=UB 

27  COHTIIUE 
12  COHTIIUE 
RETURM 
EHD 

C*************#********>*****************************SUBROUTIHE  ZLOAD. 
C  REFEREHCE:  COKPUTATIOI  OF  RADIATIOI  AID  SCATTERIIG  FOR 
C  LOADED  BODIES  OF  REVOLUTIOI 

C  HARRIIGTOH  AID  MAUTZ 

C  REPORT  AFCRL-70-0046 

C  AIRFORCE  CAMBRIDGE  RESEARCH  UBS 

C  COITRACT  10.  F-19628-68-C-0180 

SUBROUT I IE  ZLOAD ( IP, RH,ZH,ZO,Z) 

C 

C  COMPUTES  IMPEDAICE  MATRIX  ELEMEITS  FOR  LOADED  BODIES  OF  REV 
C  ZO(I)  IS  THE  SURF  IMPEDAICE  OF  THE  ITH  SEGMEIT  (IP-1  SEGMEITS) 

C  Z(.)  ARE  THE  IMPEDAICE  MATRIX  TERMS  (TRIDIAGOIAL  FOR  T-T 
C  SUBMATRIX;  DIAGOIAL  FOR  P-P  SUBMATRIX).  STORED  II  COL  VECTOR. 

C 

COMPLEX  Cl ,C2,Z0(400) , Z(2400) , XI ,X2,X3,Y1 ,Y2,Y3,FI(400) 

COMPLEX  U1,U2,U3,XI,YI 

DIMEISIOI  RH(400) ,ZH(400) ,RS(400) ,D(400) ,SV(400) 

PI=3. 14159 

MT=IP-2 

MP=IP-1 

I=MT+MP 

DO  10  IP=2,IP 

II=IP-1 

DR=RH(IP)-RH(II) 
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DZ=ZH ( IP ) -ZH ( II ) 

D(II)=SQRT(DR*DR+DZ*DZ) 

SV(II)=DR/D(II) 

RS(II)=.5*(RH(IP)+RH(II) ) 

DS=D(II)*SV(II)/2. 

Q1=RS(II)+DS 
Q2=RS(II)-DS 
FM(II)=1 . 

IF( ( ABS (Q2 ) . GT . 1 . E-6 ) .AMD. (ABS(qi) .GT. l.E-6)) 

*  FH(II)=AL0G(qi/Q2) 

10  COMTIMUE 
L0=MT*3-2 
DO  20  1=1, MP 
C1=PI*Z0(I) 

IF(I.EQ.MP)  GO  TO  80 
KI=2 

IF(I.Eq.l)  KI=1 
IF(I.Eq.MT)  KI=3 
11=1+1 

C2=PI*Z0(II) 

A=SV(I) 

IF(ABS(A) .LT. l.E-6)  GO  TO  41 
Xl=Cl*FI(I)/2./A 

X2=C1*2./A*(1.-RS(I)*FI(I)/D(I)/A) 

X3=-X2*RS(I)/D(I)/A 
GO  TO  42 

41  COITIIUE 
Xl=Cl/2./RS(I)*D(I) 

X2=(0. ,0. ) 

X3=Cl*D(I)/6./RS(I) 

42  COMTIMUE 
A=SV(II) 

IF(ABS(A).LT. l.E-6)  GO  TO  45 
Yl=C2*FH(II)/2./A 

Y2=C2*2./A*(1.-RS(II)*FM(II)/D(II)/A) 
Y3=-Y2+RS(II)/D(II)/A 
GO  TO  40 
45  COMTIMUE 

Yl=C2/2./RS(II)*D(II) 

Y2=(0 . ,0.) 

Y3=C2*D(II)/6. /RS(II) 

40  COITIIUE 
C 

C  DEFIME  TRIDIAGOMAL  ELEHEITS  FOR  T-T  SUBNATRIX  (STORED  II  COLS) 
C  (Ul-  DIAG;  U2-  LOVER;  U3-  UPPER) 

C 

XI=X1+X2+X3 

YI=Y1-Y2+Y3 

IF(RH(I).LT. l.E-6)  XI=C1/SV(I) 

IF(RH(II).LT. l.E-6)  YI=C2/SV(II) 

U1=XI+YI 
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U2=X1-X3 

U3=Y1-Y3 

L=2+(I-2)*3 

IF(KI.EQ.l)  L=0 

L1=L+1 

L2=L+2 

L3=L+3 

go  to  (50,60,70) ,ki 
50  Z(L1)=U1 
Z(L2)=U2 
GO  TO  80 
60  Z(L1)=U3 
Z(L2)=U1 
Z(L3)=U2 
GO  TO  80 
70  Z(L1)=U3 
Z(L2)=U1 

80  Z(L0+I)=2. *C1*D(I)/RS(I) 

20  CONTINUE 
RETURI 
END 

SUBROUTINE  ZTOT(MT,MP,ZL,Z) 

C 

C  ADDS  THE  SURF  IMPEDANCE  TERMS  TO  THE  TRIDIAGONAL  ELEMENTS  OF 
C  THE  BOR  IMPEDANCE  MATRIX  Z. 

C 

COMPLEX  ZL ( 2400 ) . Z ( 100000 ) 

H=MT+MP 

M0=MT*3-2 

DO  100  1=1. MP 

L0=MT*N+(I-1)*N+HT 

IF(I.Eq.MP)  GO  TO  80 

KI=2 

IF(I.Eq.l)  KI=1 

IF(I.Eq.MT)  KI=3 

L2=(I-1)*N+I 

L1=L2-1 

L3=L2+1 

M=2+3*(I-2) 

IF(KI.Eq.l)  M=0 
M1=M+1 
M2=M+2 
M3=M+3 

go  to  (50,60,70) ,ki 
50  Z(L2)=Z(L2)+ZL(M1) 

Z(L3)=Z(L3)+ZL(M2) 

GO  TO  80 

60  Z(L1)=Z(L1)+ZL(M1) 

Z(L2)=Z(L2)+ZL(M2) 

Z(L3)=Z(L3)+ZL(M3) 

GO  TO  80 
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70  Z(L1)=Z(L1)+ZL(M1) 

Z(L2)=Z(L2)+ZL(M2) 

80  Z(L0+I)=Z(L0+I)+ZL(M0+I) 

100  CONTINUE 
RETURH 
EHD 

O***********************************************SUBR0UTINE  OGIVE 
C  SUBROUTINE  :  OGIVE 
C  DATE  :  4  SEPTEMBER  1991 

C  PROGRAMMER  :  R.M.  FRANCIS 
C  REVISED  :  26  JANUARY  1992 
C  COMMENTS 

C  THIS  SUBROUTINE  WILL  GENERATE  DATA  FOR  A  BODY  OF  REVOLUTION  (BOR)  IN 
C  THE  FORM  OF  AN  OGIVE. 

C  DIMENSIONS  ARE  NORMALIZED  TO  WAVELENGTH,  SEE  PAGE  14  FRANCIS  THESIS. 

C  ZH  =  Z  CO-ORDINATE  *  2*PI 
C  RH  =  RADIUS  *2*PI 

SUBROUTINE  OGIVE ( HP , ARAD , ZH , RH , B ASE , RS , ZP ) 

C  NP  =  NUMBER  OF  POINTS  ON  THE  OGIVE  SURFACE,  MAXIMUM  =  1000 
C  ZP  =  ZPRIME,  THE  POSITION  ON  Z  WHERE  THE  RADIUS  OF  CURVATURE  STARTS 
C  BASE  »  BASE  RADIUS 

C  RS  =  RADIUS  OF  CURVATURE  IN  THE  RZ,Z  PLANE  OF  THE  OGIVE 
INTEGER  I.NP.NPBASE 

REAL  RH ( 400 ), ZH ( 400 ), ZP , BASE , RS , ZCOORD , RAD IUS , AL 
PI=3. 1415926 

C  INPUT  THE  VARIABLES  FOR  THE  OGIVE, ZP.B.RS, HP 

WRITEC6,*)  'ENTER  SURFACE  CURVATURE  (wavelengths)’ 

READC5,*)  RS 

WRITE(6,*)  ’ENTER  ZPRIME, WHERE  CURVATURE  STARTS  (wavelengths)’ 
READ ( 5 , * )  ZP 

WRITE(6,*)  ’ENTER  BASE  RADIUS  (wavelengths) ’ 

READ(5 , *)  BASE 
C  PERFORM  CALCULATIONS 

AL=SQRT(2 . *BASE*RS-BASE**2 ) 

ZMAX=  AL  +  ZP 
AHG=ASIH(AL/RS) 

L=AHIHT( (RS* AHG+ZP ) *5 ) 

HP=2*L+1 

WRITE(6,*) 'NUMBER  OF  POINTS  IS:  \NP 
DZ=  ZMAX/FL0AT(NP-1) 

DO  10  1=1, NP 

ZCOORD=ZMAX-  FL0AT(I-1)*DZ 
IF (ZCOORD. EQ. 0. ) THEN 
ZC00RD=0 . 0000000001 
ENDIF 

ZH(I)=2.*PI*ZC00RD 

IF  ( ZCOORD. LE.ZP)  THEN 
RADIUS=BASE 

ELSE 

RADIUS=SQRT(RS**2-(ZC00RD-ZP)**2)+(BASE-RS) 
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IF(RADIUS.EQ.O. )THEH 
RADIUS=0. 0000000001 
END  IF 
ESDIF 

RH ( I ) =2 . *PI*RADIUS 
10  CONTINUE 
RETURN 

:nd 
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SUBROUTINE  :  GENEX 
DATE  :  14  JANUARY  1992 

REVISED  :  17  February  1992 
PROGRAMMER  :  R.M.  FRANCIS 

SUBROUTINE  GENEX (MODE , NP , NT , NPHI , CNPHI , XT , AT , X , A , 

*  XP , AP , THS , PHS , ARAD , RBO , ZHO , R) 

SOURCE  OH  Z  AXIS  AT  Z=0.  THIS  VERSION  DOES  THE  ANTENNA  INTEGRATION 
IN  RECTANGULAR  COORDINATES.  CNPHI  POINTS  USED  IN  X  AND  Y. 
**************************** 


t  t  i 

(  V  )  =  <  W  ,  E  > 

n  i  ni 


phi  phi  i 

(  V  )  =  <  V  .  E  > 

n  i  ni 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


C* **************************** 

C  THE  EXCITATION  VECTOR  IS  COMPUTED  FOR  THE  GIVEN  R, THETA  AND  PHI 
C  COMPONENTS  ON  SURFACE  OF  BOR,  SPECIFIED  IN  SUBROUTINE  ClhCRTP. 
COMPLEX  R(800) , CEXP , PSI , SI , S2 , S3 , S4 , S5 
COMPLEX  CIRCR , CIRCT , CIRCP 

DIMENSION  RH(400) ,ZH(400) ,XT(4) , AT(4) ,X(100) ,A(100) 
DIMENSION  XP(IOO) ,AP(100) ,RH0(400) ,ZH0(400) 

INTEGER  CNPHI 
PI=3. 141592654 
BK=2 . *PI 
MP=NP-1 
MT=MP-1 
N=MT+MP 

C  LIMITS  ON  PHI  INTEGRATION 
Pl=(2.*PI-0. )/2. 

P2=(2.*PI+0. )/2. 

DO  10  1=1, NP 
RH(I)=RHO(I)/BK 
10  ZH(I)=ZHO(I)/BK 
DO  30  IP=1 ,MP 

C  QUANTITIES  FOR  THE  FIRST  SEGMENT  (POSITIVE  SLOPE) 

I=IP+1 

II=IP 
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DR=RH(I)-RH(II) 

DZ=ZH(I)-ZH< II) 

D 1 =SQRT (DR*DR+DZ*DZ ) 

Rl=. 5*(RH(I)+RH(II) ) 

Zl= . 5*(ZH(I)+ZH(II) ) 

SVP1=DR/D1 

CVP1=DZ/Di 

V1=ATAN2(DR,DZ+1 .2-5) 

C  QUANTITIES  FOR  THE  SECOND  SEGMENT  (NEGATIVE  SLOPE) 

C  (SKIP  IF  LAST  SEGMENT) 

I=IP+2 

II=IP+1 

DR=RH(I)-RH(II) 

DZ=ZH(I)-ZH(II) 

D2=SQRT (DR*DR+DZ*DZ ) 

R2=. 5*(RH(I)+RH(II) ) 

Z2=. 5*(ZH(I)+ZH(II) ) 

SVP2=DR/D2 

CVP2=DZ/D2 

V2=ATAN2(DR.DZ+1 .E-5) 

C  BEGIN  PHI  INTEGRATION:  R  TERMS  IN  Si  AND  S2;  THETA  TERM  IN  S3  AND  S4; 
C  PHI  TERM  IN  S5. 

S1=(G. ,0. ) 

S2=(0 . ,0.) 

S3=(0 . ,0.) 

S4=(0. ,0.) 

S5=(0 . ,0.) 

DO  20  J=1,NPHI 
PH=P1*X(J)+P2 

P3I=CEXP ( CMPLX ( 0 . , -MODE*PH ) ) *A ( J ) 

IF(IP.LE.MT)  THEN 

r.  t-CURRENT  INTEGRATION  FOR  THE  POSITIVE  SLOPE 
C  Gauss  quadrature 
DO  13  L=1 ,MT 
TP=n*XT(L)/2. 

RHB=R1+TP*SVP1 
ZHB=Z1+TP*CVP1 
TH=ATAN2(RHB ,ZHB+1 . E-5) 

CC=C0S(V1-TH) 

SS=SIN(V1-TH) 

CALL  CIRCRTP ( CNPHI , XP , AP , ARAD , THS , PHS , 

*  PH , RHB , ZHB , CIRCR , C IRCT , CIT  rn ) 

S1=S1+AT(L)*(.5+TP/D1)*CC*CIRCR*PSI 
S2=S2+AT(L)*( . 5+TP/D1 )*SS*CIRCT*PSI 
13  CONTINUE 

C  t-CURRENT  INTEGRATION  FOR  THE  NEGATIVE  SLOPE 
C  Gauss  quadrature 
DO  14  L=1,NT 
TP=D2*XT(L)/2. 

RHB=R2+TP*SVP2 

ZHB=Z2+TP*CVP2 
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TH=ATAI2(RHB ,ZHB+1 . E-5) 

SS=SII(V2-TH) 

CC=C0S(V2-TH) 

CALL  CIRCRTP (CIPHI , XP , AP , ARAD , THS , PHS , 

*  PH,RHB,ZHB,CIRCR,CIRCT,CIRCP) 
S3=S3+AT(L)*(.5-TP/D2)*CC*CIRCR*PSI 
S4=S4+AT(L)*(.5-TP/D2)*SS*CIRCT*PSI 

14  COITIIUE 
EHDIF 

C  phi-CURREIT  IITEGRATIOI 
DO  15  L=1 ,HT 
TP=Dl*XT(L)/2. 

RHB=R1+TP*SVP1 
ZHB=Z1+TP*CVP1 
TH=ATAI2 ( RHB , ZHB+1 . E-5 ) 

CALL  CIRCRTP ( CIPHI , XP , AP , ARAD , THS , PHS , 

*  PH , RHB , ZHB , CIRCR , CIRCT , CIRCP ) 
S5=S5+AT(L)*CIRCP*RHB/R1*PSI 

15  COITIIUE 

20  COITIIUE 

C  COMPOIEITS  ARE  STORED  II  A  COLUMI  VECTOR:  VT(1,MT),  VP(MTM.I) 
IF(IP.LE.MT)  THE! 

R(IP)=(Sl+S2)*Dl/2. *Pl+(S3+S4)*D2/2. *P1 
EIDIF 

R(IP+MT)=SS*Pl*Dl/2. 

30  COITIIUE 

RETURI 
EID 


O**  ******  *****************************************SITBR0UTIIE  CIRCRTP 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTIIE 

DATE 

REVISED 

PROGRAMMER 

2 ID  REVISIOI 

PROGRAMMER 


CIRCRTP 
1  OCTOBER  1991 

26  February  1992 
R.M.  FRAICIS 

27  October  1992 
K.  A.  KLOPP 


COMMEITS  :  THIS  SUBROUTIIE  WILL  RIGOROUSLY  CALCULATE  THE  ELECTRIC  FIELD 
FOR  A  CIRCULAR  APERTURE.  THE  FIELD  IS  CALCULATED  AT  THE  COORDIIATES 
SPECIFIED  BY  RH( . )  AID  ZH( . )  THE  APERTURE  IS  LOCATED  AT  Z  =  0, 

AID  SCAIIED  TO  A  DIRECTIOI  (THS, PHS).  THE  SUBROUTIIE  OGIVE  IS  THE 
SOURCE  OF  THE  GEOMETRIC  DATA  REQUIRED  BY  CIRC  TO  PERFORM  COMPUTATIOIS . 
ALL  PHYSICAL  DIMEISIOIS  ARE  IORMALIZED  TO  VAVELEIGTH. 


+ant  +ant 

i  \  \  rl 

E  =  -jn/(4*pi*k)  /_  /_  e  ( 

x*0  y*0 


)*TAPER(rp) 
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SUBROUTINE  CIRCRTP(CNPHI , XP , AP , ARAD , THS , PHS , 

*  PHI ,RZ,Z,CIRCR,CIRCT,CIRCP) 

INTEGER  CNPHI 
DIMENSION  XP(IOO) ,AP(100) 

REAL  K 

COMPLEX  J , JK , G1 , G2 , XI , CON , CC , SUMX , SUMY , SUMZ 
COMPLEX  CIRCR,CIRCT,CIRCP 
PI=3. 141592654 
ETA=120.*PI 

C  SET  FIELD  COMPONENTS  TO  ZERO  IN  THE  REAR  HEMISPHERE 
CIRCR=(0. ,0.) 

CIRCT=(0. ,0. ) 

CIRCP=(0. ,0. ) 

C  Z  OF  THE  OBSERVATION  POINT  IS  >=  0  IN  THE  FORWARD  HEMISPHERE 
IF(Z.GE.O. )  THEN 
SUMX=(0. 0,0.0) 

SUMY=(0. 0,0.0) 

SUMZ=(0. 0,0.0) 

K=6. 283185307 
J  =  (0.0, 1.0) 

JK=  (0.0,6.283185307) 

C0N=-J/(4.*PI) 

C  OMIT  THE  FACTOR  OF  ETA  SINCE  IT  IS  NOT  INCLUDED  IN  SUBROUTINE  ZMAT. 
C  ARAD**2  IS  FOR  SCALING  OF  GAUSSIAN  INTEGRATION.  READJUST  SCALING 
C  TO  AGREE  WITH  OLD  SUBROUTINE  SPHERE  (ARBITRARY) 

CC=  CON 

C  OUTER  LOOP:  INTEGRATE  OVER  X,  -ARAD  <  X  <  ARAD. 

C  INNER  LOOP:  INTEGRATE  OVER  Y,  -ARAD  <  Y  <  ARAD. 

C  EVALUATE  ANTENNA  FIELD  AT  RZ,Z: 

XO=RZ*COS(PHI) 

YQ=RZ*SIN(PHI) 

STS*SIN(-THS) 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 

C  INTEGRATE  IN  X.Y  COORDINATES 
DO  50  M=l, CNPHI 
X=ARAD*XP(M) 

DO  60  N*l, CNPHI 
Y=ARAD*XP(N) 

RP=SQRT(X**2+Y**2) 

C  SKIP  INTEGRATION  POINTS  BEYOND  THE  ANTENNA  RADIUS 
IF(RP.LE.ARAD)  THEN 
PHIPRIME=ATAN2 ( Y , X+ 1 . E- 10 ) 

R=SqRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 
G1=(((K*R)**2)-1-(JK*R))/R**3 
G2=(3+(3*JK*R)-(K*R)**2)/R**S 
X 1 = JK* (RP*COS(PHIPRIME-PHS ) *STS-R) 

Y1=X0-X 

X2=Y1**2 

Y2*Y0-Y 

SUMX»SUMX+(CC*AP(M)*AP(N)*(G1+(X2*G2))*CEXP(X1))*TAPER(RP) 
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SUMY=SUMY+(CC*AP(M)*AP(N)*Y1*Y2*G2*CEXP(X1) )*TAPER(RP) 
SUMZ=SUMZ+(CC*AP(M)*AP(N)*Y1*Z*G2*CEXP(X1))*TAPER(RP) 
ENDIF 

60  CONTINUE 
50  CONTINUE 

C  ASSUME  AI  ELEMENT  FACTOR,  ELMT 
c  ELMT=SqRT(ABS(C) ) 

ELMT=1 . 


CIRCT=(C*COS(PHI)*SUMX+C*SIN(PHI)*SUMY-S*SUMZ)*ELMT 

CIRCR=(S*CQS(PHI)*SUMX+S*SIN(PHI)*SUMY+C*SUMZ)*ELMT 

CIRCP= ( -SII (PHI ) *SUMX+COS (PHI ) *SUMY ) *ELMT 

EIDIF 

RETURI 

END 


C**********************************************SUBR0UTINE  TESTSPHERE 
C  SUBROUTIME  TESTSPHERE 

C  DATE  21  FEBRUARY  1992 

C  PROGRAMMER  R.  M.  FRANCIS 

C  COMMENTS  GENERATES  A  SPHERE,  WITH  PLOTTED  POINTS 

C  AT  EQUAL  INTERVALS  IN  THETA. 

SUBROUTINE  TESTSPHERE(NP , ZH , RH , BASE , RS , ZP ) 

INTEGER  NP,I 

REAL  ZH(400) ,RH(400) . BASE , RS , ZP , SPRAD 
PI=3. 1415926 


BK=2.*PI 

ZP=0.0 

WRITE (6,*) 'ENTER  NUMBER  OF  POINTS  (IP):’ 
READ(5,*)NP 

WRITE (6 , * ) ' ENTER  SPHERE  RADIUS ' 

READ (5,*) SPRAD 
BASE=SPRAD 
RS=SPRAD 
DO  1241  1=1, NP 

AIGLE=PI*FLOAT( 1-1 ) / ( 2 . * FLOAT ( I P- 1 ) ) 
ZH ( I ) =BK*SPRAD*COS ( ANGLE ) 

RH ( I ) =BK*SPRAD*SIN (ANGLE ) 

IF(RH(I) .Eq.O. )THEN 
RH(I >=0.0000000001 
ENDIF 

1241  CONTINUE 
RETURN 
END 


C«*******«******************«*****#«**«**«*«*«**«**SUBROUTIIE  CONE 
C  GENERATES  THE  BOR  GEOMETRY  FOR  A  CONE. 

SUBROUTINE  COIE(IP,ZH,RH,BASE,RS,ZP) 

REAL  ZH(400) , RH(400) ,BASE,RS,ZP 

REAL  HA , AIG , ZMAX , DZH 

INTEGER  NP 

PI=3. 1415926 

BK=2 . *PI 

WRITE (6,*)’ ENTER  THE  CONE  HALF  ANGLE  IN  DEGREES' 
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READ(5 , *)HA 

IF(HA . GT. 90 . )THEI 

WRITE ( 6 ,  * )  ’  RAIGE  OF  AHGLE:  0<  AHGLEOO  ’ 

WRITE ( 6 , * ) 1 RAIGE  EXCEEDED  * 

GOTO  10 
EIDIF 

ANG=HA*PI/180. 

WRITE ( 6 , * ) ’ EITER  THE  MAXIMUM  Z’ 

READ(5,*)ZMAX 

WRITE (6,*) ’EITER  THE  lUMBER  OF  POUTS  (IP)’ 

READ(5 , *)IP 
ZH(1)=ZMAX*BK 
RH ( 1 ) =0 . 

DZH=ZH ( 1 ) /FLOAT ( IP - 1 . ) 

DO  10  1=1, IP 

IF(ZH ( I ) . EQ . 0 . ) THE! 

ZH (I) =.00000001 
EIDIF 

ZH(I+1)=ZH(I)-DZH 

RH(I+1)=RH(I)+(DZH*SII(AIG)) 

10  COITIIUE 

BASE=RH(IP)/BK 

ZP=0. 

RS=0. 

RETURI 

EID 

O****’M‘=********************************’M*****SUBR0UTIIE  DISK 
SUBROUTIIE  D I SK ( IP , ZH , RH , B ASE , RS , ZP ) 

REAL  ZH(400) , RH(400) , BASE , RS , ZP , RMAX 
IHTEGER  IP 

WRITE (6,*) ’EITER  THE  RADIUS  OF  THE  DISK’ 

READ (5,*) RMAX 

WRITE (6,*) ’EITER  THE  DISTAICE  TO  THE  DISK’ 

READ(5,*)ZP 

WRITE (6,*) ’EITER  THE  IUMBER  OF  POUTS  (IP)’ 

READ(5,*)IP 
PI=3. 1415926 
BK=2 . *PI 

DRH=RMAX*BK/FL0AT(IP-1. ) 

DO  10  1=1, IP 
ZH(I)=ZP*BK 
RH ( I ) =DRH*FLOAT ( I- 1 . ) 

10  COITIIUE 
BASE=RMAX 
RS=0 . 

RETURI 

EID 

O******«*********************************«****SUBR0UTIIE  PARAB 
SUBROUTIIE  PARAB ( IP , ZH . RH , DM , FOD , ZP ) 

DIMEHSIOI  RH(400) ,ZH(400) 

WRITE (6,*)  'EITER  DIAMETER  AID  F/D  RATIO’ 
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READ(5,*)  DM.FOD 

WRITE(6,0  ’FEED  DISTAICE  FROM  FOCUS  (+  IS  VEARER) ’ 
READ(5,*)  ZP 

WRITE (6.*)  ’ESTER  HUMBER  OF  GEIERATIIG  POUTS’ 

READ(5 , *)  IP 
BK=2 . *3 . 14159 
FM=FQD*DM 

PHIV=2.*ATAH(l./4./F0D) 

DO  24  1=1, HP 

TH=FLOAT( I- 1 ) ♦PHIV/FLOATdP- 1 ) 

RM=2 . *FM/ ( 1 . +COS (TH ) ) ♦BK 
ZH  ( I )  =RM*COS  (TH )  +ZP*BK 
RH(I)=RM*SII(TH) 

24  COITIHUE 
RETURH 
EHD 

^♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦FUICTIOI  TAPER 
FUICTIO*  TAPER(RHO) 

C  SPECIFY  AHTEHIA  ILLUMII-TIOI  FUICTIOH.  REAL  FUSCTIOI  OF 
C  RHO  OILY 

TAPER=1 . 

RETURH 

EHD 

(^♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦SUBROUTHE  AITFF 

SUBROUTIIE  AITFF (THETA , PHI , THETAS , PHIS , ARAD , ETF , EPF ) 

C  COMPUTES  CLOSED  FORM  AITEIIA  PATTER!  II  THE  FAR  FIELD. 

C  ETF  AID  EPF  ARE  RETURIED.  AIGLES  ARE  II  RADIAIS. 

C  PRESEHT  EITRY  IS  FOR: 

C 

C  UHIFORM  ILLUMIHATIOI  SCAIIED  TO  A  DIRECTION  (THETAS, PHIS) 
C 

COMPLEX  ETF , EPF , JK , SCL , EB 
JK=(0. 0,6. 283185307) 

BK=6. 283186307 
PI=BK/2 . 

C  SET  RADIATIOI  PATTER*  TO  ZERO  I*  THE  REAR  HEMISPHERE 
IF(COS (THETA) .LT.O. )  THE! 

EB=(0. ,0.) 

ELSE 

BB=SII (THETA) ♦SII(PHI )-SII (THETAS ) *SII (PHIS ) 

AA=SII (THETA ) *COS(PHI )-SII (THETAS )*COS( PHIS) 
ARG=BK*ARAD*SQRT(AA*^2+BB**2) 

SCL=-(0. , 1 . )*2. ♦PI**2 
IF(ABS(ARG) .LT. 1 .E-5)THEI 
EB=SCL/2. 

ELSE 

EB=SCL* (BESS J 1 ( ARG ) / ARG ) 

EIDIF 

EIDIF 

C  EB  IS  THE  X  COMPOIEIT ;  GET  ETHETA  (ETF)  AID  EPHI  (EPF) 

C  ASSUME  A!  ELEMEIT  FACTOR,  ELMT 
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ELMT=SQRT(ABS(COS (THETA) ) ) 

ELMT=1. 

ETF=  EB  *COS ( THETA ) *COS ( PHI ) *ELMT 
EPF=-EB*SIH(PHI)*ELMT 
RETURN 
END 

FUNCTIOH  BESSJ(N.X) 

C  RETURNS  THE  BESSEL  FUHCTION  B  OF  ORDER  N  (>l)  AND  REAL 
C  ARGUMENT  X. 

PARAMETER  (IACC=40 , BIGNO=l . ElO , BIGNI=1 . E-10) 
IF(N.EQ.O)  THEN 
C  IF  N=0  CALL  BESSJO 

BESSJ=BESSJO(X) 

ELSE  IF(N.Eq.l)  THEN 
C  IF  N=1  CALL  BESSJ1 

BESSJ=BESSJ1(X) 

ELSE 

C  IF  N>1  USE  RECURSION 
BESSJ=0 . 

IF(X.NE.O.)  THEN 
T0X=2./X 

IF(X.GT.FLOAT(N) )  THEN 
BJM=BESSJO(X) 

BJ=BESSJ1(X) 

DO  11  J=1,N-1 

BJP=J*TOX*BJ-BJM 

BJM®BJ 

BJ=BJP 

11  CONTINUE 
BESSJ=BJ 
ELSE 

M=2*( (N+INT(SQRT(FLOAT(IACC*N) ) ) )/2) 

8ESSJ=0. 

JSUM=0 . 

SUM=0 . 

BJP=0. 

BJ=1 . 

DO  12  J=M ,1,-1 

BJMsJ*TOX*BJ-BJP 

BJP=BJ 

BJ=BJM 

IF(ABS(BJ) . GT.BIGNO)  THEN 

BJ=BJ*BIGNI 

BJP*BJP*BIGNI 

BESS J=BESS J*BIGNI 

SUM=SUM*BIGNI 

ENDIF 

IF(JSUM.NE.O)  SUMaSUM+BJ 
JSUM=1-JSUM 
IF(J.EQ.N)  BESS J=BJP 

12  CONTINUE 
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SUM=2 . *SUM-BJ 

BESSJ=BESSJ/SUM 

ENDIF 

ENDIF 

END  IF 

RETURN 

END 

FUNCTION  BESSJO(X) 

C 

C  BESSEL  FUNCTION  OF  0  ORDER,  REAL  ARGUMENT  X 
C  (SEE  ’NUMERICAL  RECIPES’,  P.172) 

C 

REAL*8  Y,P1 ,P2,P3 ,P4,P5,Q1 ,Q2 , Q3 ,Q4, Q5 ,R1 ,R2,R3,R4,R5,R6 , 

*  SI ,S2 ,S3 ,S4 ,S5 ,S6 

DATA  PI ,P2,P3 ,P4,P5/1 .DO , - . 109862827D-2 , . 27345 10407D-4 , 

*  - . 2073370639D-5 , . 20938872 1 ID-6/ 

DATA  Ql ,Q2 ,Q3,Q4,Q5/~ . 156249999SD- 1 , . 1430488765D-3 , 

*  -.6911 14765 ID-5 , . 7621095161D-6 , - . 934945152D-7/ 

DATA  R1 , R2 , R3 , R4 , RS , R6/S7568490574 . DO ,-13362590354 . DO , 

*  651619640 . 7D0 , -11214424 . 18D0 , 77392 . 33017D0 , -184 . 9052456D0/ 
DATA  SI , S2 , S3 , S4 , SS , S6/5756849041 1 . DO , 1029532985 . DO , 

*  9494680 . 7 18D0 , 59272 . 64853D0 , 267 . 8532712D0 , 1 . DO/ 

IF(X.Eq.O.)  THEN 

BESSJ0=1. 

ELSE  IF(ABS(X) .LT.8. )  THEN 
Y=X**2 

BESS  J0=  (R1+ Y*(R2+Y*  (R3+Y*  (R4+Y*  (R5+Y*R6  )))))/ 

*  (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ) ) ) 

ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX-. 785398164 

BESSJO=SQRT( . 636619772/AX)*(C0S(XX)*(P1+Y*(P2+Y*(P3+ 

*  Y*(P4+Y*P6) ) ) ) -Z*SIN(XX) * (qi+Y*(q2+Y* (Q3+ 

*  Y*(q4+Y*q5))))) 

ENDIF 

RETURN 

END 

FUNCTION  BESSJl(X) 

C 

C  BESSEL  FUNCTION  B  OF  ORDER  1,  REAL  ARGUMENT  X 
C  (SEE  ’NUMERICAL  RECIPES',  P.173) 

C 

REAL*8  Y,Pl,P2,P3,P4,P5,qi,q2,q3,q4,qS,Rl,R2,R3,R4,R5,R6, 

*  SI ,S2,S3,S4,S5,S6 

DATA  P1,P2,P3,P4,P5/1.D0, . 183105D-2 , - . 35 16396496D-4 , 

*  .24S7520174D-6, -.203370190-6/ 

DATA  qi , q2 , q3 , q4 , qs/ . 04687499995D0 , - . 2002690873D-3 , 

*  . 8449 199096D-5 , - . 99228987D-6 , . 1057874 1 2D-6/ 

DATA  R1 ,R2,R3,R4,R5, R6/723626 14232 . DO , -7895059235 . DO , 
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*  242396853 . IDO , -29726 1 1 . 439D0 , 15704 . 4826D0 ,-30.1 6036606D0/ 
DATA  51,32,33,34,35, S6/144725228442 . DO , 2300535178 . DO , 

*  18583304 . 74D0 , 99447 . 43394D0 , 376 . 9991397D0 , 1 . DO/ 
IF(X.EQ.O.)  THEH 

BESS J 1=0. 

ELSE  IF(ABS(X) .LT.8. )  THEH 
Y=X**2 

BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(RS+Y*R6)))))/ 

*  (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ) ) ) 

ELSE 

AX=ABS(X) 

2=8. /AX 
Y=Z**2 

XX=AX-2. 356194491 

BESSJ1=SQRT( . 636619772/AX)* (COS (XX)* (P1+Y*(P2+Y*(P3+ 

*  Y*(P4+Y*P5) ) ) )-Z*SIH(XX)*(qi+Y*(Q2+Y*(q34- 

*  Y*(q4+Y*qS)))))*SIGI(l. ,X) 

EHDIF 

RETURI 

EID 
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APPENDIX  C. 


GAIN.F  LISTING 


Material  on  the  following  pages  is  a  listing  of  the  program  GAIN.F  described 
in  chapter  III. 
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c******************** ************************************ ••*••*•»•***••* 

c 

C  PROGRAM  :  gain.f 

C  DATE  :  17  November  1992 

C  PROGRAMMERS  :  Dr.  D.  Jenn,  LT  K.  Klopp 
C 

C  Computes  the  gain  of  a  circular  aperature  antenna  radiating  through 
C  a  radome  consisting  of  a  body  of  revolution,  by  numerically 
C  integrating  electric  field  surrounding  system.  E-Field  is  generated  by 
C  currents  distributed  on  a  surface  of  the  body  of  revolution.  The 
C  currents,  and  their  locations  sere  previously  determined  using  the 
C  program  "radome. f",  and  stored  in  the  file  "f77curcoef sA  or  B." 

C 

C  INPUT  <  f77curcoefsA/B 
C 

C***************************************************************»******* 
CHARACTER  ITH 
CHARACTER* 12  CC 
CHARACTER* 10  GO 
CHARACTER*8  GNT.ONPHI.CGP 
CHARACTER* 14  TPTS ,PPTS .PHIPTS 
COMPLEX  R( 1 600 ) , LEADTERM , ETT , EPT 
COMPLEX  EP , ET , ETF , EPF , ETCOMB , EPCOMB 
COMPLEX  EXP 1 , EXP2 , CON JG , CEXP , CMPLX , CT ( 30000 ) , IMP , JK 
COMPLEX  U , UC , ET 1 , EP 1 , ET2 . EP2 , ZLO (400 ) 

C  DIMENSION  RH(400) ,ZH(400) , IPS(800) 

DIMENSION  RH(400) ,ZH(400) 

DIMENSION  XT(4) ,AT(4) ,A(100) ,X(100) ,XP(100) ,AP(100) 

DIMENSION  Q(30) ,S(30) 

DIMENSION  THRR( 1000, 1000), PHRR(IOOO) 

C  DIMENSION  ETSCATC500) .EPSCAT(SOO) 

C  DIMENSION  EXP(500) , ANG(500) ,ECP(SOO) 

INTEGER  CNPHI, SELECTION 
DATA  PI/3.1415926/ 

DATA  RR/ 1000/ 

Rad=PI/180 . 

BK=2 . *PI 
U=(0. ,1.) 

U0=(0. ,0. ) 

UC=-U/4./PI 
JK=(0. 0,6. 283185307) 

C*****ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS***************** 
C*****  (ALL  DATA  FILES  ARE  APPENDED  WITH  THIS  LETTER)  *********** 

C*****  (.m  FILES  ARE  FOR  GENERATING  PLOTS  IN  MATLAB)  *********** 

WRITE(6,*)’ ENTER  A  LETTER  TO  INDICATE  WHICH  ITERATION  THIS  IS’ 

RE AD (5,*) ITH 
CC=’curcoefsdat '//ITH 
GO=’gainout ’//ITH//’ .m’ 

C********** **************************** ********************************** 

c 

C  READ  in  the  following  data  from  file  called  f77curcoefs: 
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c 

C  iteration,  ITH 

C  geometry  of  BOR,  SELECTIOI 

C  radius  of  base  of  BOR,  BASE 

C  total  points,  IP 

C  number  of  t-intervals,  MT 

C  phi-intervals,  MP 

C  sum  of  t  and  phi  intervals,  I 

C  number  of  rows  of  currents,  IROV 

C  total  number  of  modes ,  MOOES 

C  number  of  sub  vectors,  IBLOCK 

C  antenna  seem  angle  in  theta,  TBS 

C  antenna  scan  angle  in  phi,  PHS 

C  antenna  radius,  ARAD 

C  observation  angle,  PHIO 

C  complex  radome  impedance,  IMP 

C  zprime ,  ZP 

C  surface  radius,  RS 

C  closed  form  antenna  pattern,  IFF=0 

C 


C  files  share  Gauss  integration  weights  and  abscissas  are  stored 

C 

C  coordinates  of  surface  points  *2PI/lamda,  ZH(I),  and  RH 0 ) 

C 

C************************************************************************* 

0PEI(1 ,file=CC) 

WRITE (6 ,*)  'FILE  OPEIED  IS  ' ,CC 
READ (1 , • )  ITH , SELECTIOI , BASE , IP , MT , MP , I 
READ ( 1 ,  * )  IROV , MODES , IBLOCX , THS , PHS , ARAD , PHIO 
READ(1,*)  IMP , ZP , RS , IFF 
c  writ e(6,*) IMP, ZP.RS, IFF 

READ(1,*)  (RH(L),L=1,IP) 
c  write(6,*)  (RH(L) ,L=1 .IP) 

READ(1,*)  (ZH(L) ,L=1 ,IP) 
c  write(6,*)  (ZH(L) ,L=1 ,IP) 

READ(1,*)  (ZLO(L) ,L=i ,IP) 

C  write(6,*)  (ZLO(L) ,L=1,IP) 

C* . . 

c 

C  READ  in  current  coefficients,  CT(I)  from  file  curcoefsdatA/B 
C  (This  is  one  one  collumn  vector.) 

C 

CM»***««MM«M*M**MM*M*MMM*«MmM*MtM*M*****M*M«****M**» 

READd,*)  (CT(L)  ,L“1 ,  IROV) 

C  write(6,*) (CT(L) ,L=1 ,IR0W) 

READ(1,«)  GIT.GIPHI , CGP 

WRITE (6,*)  ’RUI  CODE  READ  FROM  DISC  :  ’ ,ITH 

THETAS=THS*RAD 

PHIS=PHS*RAD 

MHI=M0DES+1 

CLOSE (1) 
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c********************************************************************* *«***•*•• 

c 

C  Read  in  weights  A(I),  and  Abscissas  X(I)  for  the  following  Gauss- 
C  Quadrature  Integrations: 

C  file  wts  absc 

C 

C  BOR  integration  in  t  direction:  TPTS  =  gaus/gaus#;  XT(I),AT(I) 

C  BOR  integration  in  phi  direction:  HPHI  =  gaus/gaus#;  X  (I), A  (I) 

C  Antenna  integration  in  x,y  direction:  PHIPTS  =  gaus/gaus#;  XP(I),AP(I) 

C 

C****** **************** *********** *#*********»»<■  ************ •**•****•****•••••* 

TPTS= ’ gaus/ ’ //GST 

PPTS= ’ gaus/ ’ // GSPHI 

PHIPTS= ’ gaus/ ’ // CGP 

OPES ( 2 , FILE=TPTS , STATUS  = ' OLD ■ ) 

OPES ( 3 , FILE=PPTS , STATUS  = ’ OLD  1 ) 

OPES ( 8 , FILE=PHIPTS , STATUS= 1 OLD 1 ) 

READ ( 2 , * ) ST 
WRITE (6,*) ’ST  =  1 ,BT 
IF(ST.GT.4)THES 

WRITE (6,*) ’MAXIMUM  SUMBER  OF  POISTS(ST)  IS  4’ 

GOTO  999 
ESDIF 

READ(3 , *)BPHI 
WRITE(6,*)’SPHI  =  ’ ,SPHI 
IF(BPHI.GT.200)THES 

WRITE(6,*) ’MAXIMUM  SUMBER  OF  POIBTS(BPHI)  IS  200’ 

GOTO  999 
ESDIF 

READ (8,*) CSPHI 
WRITE(6,*) ’CSPHI  =  ’ , CSPHI 
IF(CBPHI.GT. 100)THEB 

WRITE (6.*) ’MAXIMUM  SUMBER  OF  POIBTS(CBPHI)  IS  100’ 

GOTO  999 
ESDIF 

0  LOAD  THE  WEIGHTS  ASD  ABSCISSAS  IS  THE  VECTORS. 

DO  1  M=1 ,ST 

READ(2,*,EHD=1)XT(M) , AT(M) 

WRITE(6,*)XT(M) ,AT(M) 

1  COSTISUE 

DO  2  M=1,SPHI 

READ (3 , * ,ESD=2)X(M) ,A(M) 

WRITE(6,*)X(M) ,A(M) 

2  COSTISUE 

DO  4  M*l, CSPHI 

READ(8 ,* ,ESD=4)XP(M) ,AP(M) 

WRITE(6 ,*)XP(M) ,AP(M) 

4  COSTISUE 
CL0SE(2) 

CL0SE(3) 

CL0SE(8) 
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C***********************e*********** ********* ********************************** 

c 

C  Create  an  output  file  called  ’gainoutA/B.m’ ,and  output  data  on  geometry  etc. 
C 

C****************************************************************************** 

OPEH ( 7 , FILE=GO ) 

WRITE(7 ,2000) 1  GEOMETRY  =  ' .SELECTIOH , ’  HP  =  ’.HP 


A, 

’  MT  = 

,MT, ’  KP  =  ’.MP, ’  H  = 

M 

A, 

’  MODES 

=  ’.MODES,’  IB LOCK  = 

,  IBLOCK 

A, 

’  PHI  = 

■  ,PHI 

WRITE(7,2001) ’  BASE  =  ’.BASE,’  ARAD  =  ’.ARAD,’  THETASCAH  =  ’ ,THS 
A,  ’  PHISCAH  =  ’ ,PHS 

WRITE (7,*) ’IPHI  = ’ ,HPHI 

2000  FORMAT ( / /SX ,  RADOME  GAII  CALCULATIOI  */.•/.’//, 10(/’*/.  ’.A13.I4)) 

2001  FORMAT (3 (/”/.  ’ .A13.F9.2) ,/) 

C*** ***************************************************************** *********** 

c 

C  Gain  Computation: 

C  4  pi  *  Max  Radiation  Intensity 

C  G  = _ 

C  Integration  Over  Closed  Surface  (  Radiation 

C  Intensity  *  Beam  Solid  Angle  ) 

C 

C  Radiation  Intensity  =  (l/(2*eta))*l  E(t,phi,r)  1**2 
C 

C  Since  E  is  in  the  fax  field,  need  to  only  determine  Etheta,  Ephi. 

C  _ 

C  lEl-‘*2  =  IE  scatteredl**2  +  IE  antennal**2 
C 

C  IE  scattered  I  *  Etheta  +  Ephi 
C 

C* ****************************************************************** ************ 
C**Enter  t  divisions  for  integration  in  theta,  and  phi 

C  directions.  (This  allows  easy  variation  of  the  number  of  integration  steps 
C  without  having  to  recompute  a  different  set  of  weights  and  abscissas  for 
C  Gauss  Integration. ) 

WRITE (6,*) 'EITER  HUMBER  OF  DIVISIOIS  II  THETA  DIRECTIOI’, 

A  >  IDIVT  =  ?* 

READ ( 5 , * ) IDI VT 

WRITE (6,*) ’IDIVT  =  ’.IDIVT 

WRITE(7,*)  ”/.  IDIVT  =  ’.IDIVT 

WRITE (8,*) ’EITER  HUMBER  OF  DIVISIOIS  II  PHI  DIRECTIOI,  IDIVP  =  ?’ 

READ(5 , *)IDIVP 
WRITEC6,*) 'IDIVP  =  ’ , IDIVP 
WRITE(7,*)”/.  IDIVP  =  ’.IDIVP 
kX=IDIVT 

C  do  500  IDIVT=1 ,Uc 
IDIVP=IDIVT 

C***********************************integration  intervals  0<theta<180,  0<phi<360 
THETA 1*0 
THETA2=180 
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S1=THETA1 

S2=THETA2 

Sl=Sl*Rad 

S2=S2*Rad 

PHI1=0 

PHI2=360 

Q1=PHI1 

Q2=PHI2 

qi=qi*Rad 

Q2=Q2*Rad 

C**************************»********theta  integration  angles 
DS-(S2-S1)/FL0AT(IDIVT) 

DO  200  1=1 .IDIVT+l 
S(I)=(I-1)*DS 

200  COHTIIUE 

C***********************************phi  integration  angles 
DQ=(Q2-Q1)/FL0AT(IDIVP) 

DO  201  1=1 .IDIVP+l 
q(I)=(I-l)*DQ 

201  COKTIIUE 
SUH=0. 

SUN10IE=0 . 

UMAX=0. 

UMAXB0IE=0 . 

EHAX=0. 

EMAXI0IE=0 . 

C************************V  phi  integration  of  rad  intensity  *  beast  solid  angle  V 
DO  300  JJ=1 ,IDIVP 
Pl=Dq/2. 

P2=(q(JJ+i)+q(JJ))/2. 

C**********************************V  phi  integration  determining  Escattered  V 
C  V  within  each  interval  o f  phi,  IDIVP  V 

DO  300  J=1 ,IPHI 
PHR=P1*X(J)+P2 

C *****************v  theta  integration  of  rad  intensity  *  bean  solid  angle  V 
DO  301  11=1 .IDIVT 
Tl=DS/2. 

T2=(S(II+l)+S(II))/2. 


0***«*********************v  theta  integration  determining  Escattered  V 
C  V  within  each  interval  of  theta,  IDIVT  V 


DO  301  1=1 ,IPHI 
THR=T1*X(I)+T2 
IF(THR.LT.PI)  GO  TO  302 
THR=(BK-TIR) 

PHR=PHR+PI 

302  COITIVUE 

C***************************************************Beam  Solid  Angle 
STRA=SII(THR )*A(I)*A(J) 

KJ=(JJ-1)*IPHI+J 
KI=(II-1)*1PHI+I 
THRR(KJ,KI)=THR* 180/PI 
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PHRR(KJ)=PHR*180/PI 

SIP=SI*(PHR) 

CIP=COS(PHR) 

C*****»»* ********************************************************************** 

C  _  _ 
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CV********************************see  page  25.  H,E,  and  Combined  Field********V 

C  -jkr 

C  -jhe 

C  Leading  Term  _  :  eta  (h)  is  cancelled  by  its  inverse  in  CT  (I), 

C  4pir  r  dependence  is  ignored,  and  the  shole  Term  is 

C  scaled  by  2pi  for  gauss  quadrature  integration 

C 

C  Leading  Term  is  reduced  to  (-j/4pi)/2pi 

c*************** ******************************* *******************************V 

C  LEADTERM=UC/BK 

LEADTERM=UC 
ET=(0,0) 

EP=(0 ,0) 

ET1=(0. ,0. ) 

EP1=(0 . ,0.) 

ET2=(0. ,0. ) 

EP2=(0. ,0.) 

C******************************************************************* ********** V 
C  MHI 

C 

C  \ 

C  /__ 

C  n=l 

c********* *************************************************************** *****V 

DQ  305  M*1 ,MBI 
IM=H-1 

C***** *************************************** *****************  jnphi 
C  e 

EXP1=CEXP ( CMPLX ( 0 . , FLOAT ( IH) *PHR) ) 

EXP2=C0IJG(EXP1) 

C*************determin#  plane  save  excitation  vector  in  far  field,  R 
CALL  PLA1E ( IN , IM , IP , IT , RH , ZH , XT , AT , TER , R ) 

ITOP 1 =M0DES -EM 

ITOP  2 = VBLOCX - ( ITOP 1  +  1) 

IS2*IT0P1*I 

■si=rroP2*i 

c******************************************************************* 
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c  I  _  t , theta  I  I  _  t  I  jnphi  I  _  t.phi  I  I  _  t  I  jnphi 

C  I  R  I . I  I  I  e  , and  I  R  ! . I  I  I  e 

C  In  llnl  In  llnl 

O************* ***************************************************** 

DO  306  L= 1 , MT 

C***********************************************positive  mode  +n 

ET1=ET1+R(L)*CT(L+HS1)*EXP1 

EP1=EP1+R(L+H)*CT(L+IS1)*EXP1 

IF(IM.EQ.O)  GO  TO  306 

ET1=ET1+R(L)*CT(L+IS2)*EXP2 

EP1=EP1-R(L+H)*CT(L+IS2)*EXP2 

306  COITIIUE 

C  I  _  phi.thetal  I  _  t  I  jnphi  I  _  phi, phi  I  I  _  t  I  jnphi 

C  I  R  I . I  I  I  e  , and  I  R  I . I  I  I  e 

Cln  llnl  In  llnl 

C****** ************************************************************* 

DO  307  L=1 ,MP 

ET2*ET2+R (L+MT ) *CT ( L+IS 1 +MT ) *EXP1 
EP2=EP2+R (L+MT+I ) *CT (L+IS 1 +MT) *EXP 1 
0**********************************************negative  mode  -n 
IF(IM.Eq.O)  GO  TO  307 
ET2=ET2-R (L+MT ) *CT (L+IS2+MT ) *EXP2 
EP2=EP2+R( L+MT+I ) *CT(L+IS2+HT) *EXP2 

307  COITIIUE 

c  ET=ET+LEADTERM*(ET1+ET2) 

c  EP=EP+LEADTERM* ( EP 1+EP2 ) 

305  COITIIUE 

C***********************************************'  end  mode  summation 
ETT=LEADTERH* (ET1+ET2) 

EPT=LEADTERM*(EP1+EP2) 

******************************************add  in  antenna  contribution 
C  FEED  COITRIBUTIOI  II  THE  FAR  FIELD  IS  ETF.EPF 
O **************************** ***********************BESSJ1 . 

C  USE  CLOSED  FORM  FEED  EXPRESSIOI  FROM  SUBROUTIIE  AITFF  IF 
C  IFF=0;  OTHERWISE  USE  BRUTE  FORCE  EVALUATION  FROM  CIRCRTP 


R0BS=1000. 

IF(IFF.Eq.O)  THEI 

CALL  AITFF (THR , PHR , THS , PHS , ARAD .ETF.EPF) 

C  urite(6,*)  ’  ET  =  ’ .CABS(ET) , ’  ETF  =  ’ .CABS(ETF) 

ELSE 

RHBsROBS*SII (THR) 

ZHB=R08S*C0S (THR ) 

CALL  CIRCRTP (CIPHI.XP.AP, ARAD, THS, PHS, 
k  PHR , RHB , ZHB , CIRCR , CIRCT , CIRCP ) 

C  REMOVE  THE  1/R  DEPEIDEICE  BECAUSE  EXP(-jkR)/R  IS  OMITTED  II 
C  THE  SCATTERED  FIELDS  ET  AID  EP 
ETF* CIRCT* ROBS 
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EPF=CIRCP*ROBS 

C  write(6,*) ’  ET  =  ’ .CABS(ET) , 1  ETF(ROBS)  =  ’  .CABS(ETF) 

EIDIF 

C*************************************«********total  E  theta  and  E  phi 

C  srite(6,*) ’THR( ’ , j j , j ,ii,i,IM, ’ )  =  ’,180*THR/pi 

ETCOMB=ETT+ETF 
EPCOMB=EPT+EPF 

EMAG=SQRT ( CABS ( ETCOHB ) **  2+CABS ( EPCOMB ) **2 ) 

EMAGIOIE=SQRT(CABS (ETF) **2+CABS (EPF ) **2 ) 

IF (EMAG . GT . EMAX ) THEM 
EM  AX = EM  AG 
EHAXIOIE=EMAGIO*E 
WRITE(6,*) 1  EMAG ( ’ , J , I , ’ )  =  ’.EMAG 
EID  IF 

C*********************e******************determine  max  field  intensity 
UMAG=EMAG**2 
UMAGI0IE=EMAGI0«E**2 
IF (UMAG . GT . UMAX ) THE! 

UMAX=UMAG 

THETAMAX= 1 80*THR/P I 
PHIMAX=180*PHR/PI 

C  WRITE(6, *) 1  UMAX  =  ' .UMAX 

C  WRITE(6,*)’  THETAMAX  =  ’.THETAMAX 

C  WRITE(8,*) ’  PHIMAX  =  ’.PHIMAX 

EID  IF 

IF (UMAGIOIE . GT . UMAXIOIE ) THE! 

UMAXVOIE=UMAGIOIE 
EID  IF 

C********«******«*******nultiply  integration  by  beam  solid  angle  (STRA) 
SUM=SUM+UMAG*STRA 
SUMSOSE=SUMIOIE+UMAGSOIE*STRA 
301  COITIIUE 

C*****e**********************************e*e******‘  end  theta  integration  ‘ 

300  COITIIUE 

C********************************************************'  end  phi  integration  * 
C*ee****e*e*e**e*e*******e*scale  each  summation  by  (b-a)/2  for  Gauss  integration 
PRAD=T1*P1*SUM 
PRADI0IE=T1*P1*SUMI0IE 

WRITE(6,*) ’PRAD  =  ’.PRAD 

WRITE(8,*) ’PRAD  10  RADOME  =  ’.PRADIOIE 

C*************************e*******e********************************************* 
C  End  Integration  over  closed  surface  radiation  intensity  *  beam  solid  angle 
C**********«******************************e******e*************************«**»* 
C*******************************************************************compute  gain 
GAII=4*PI*UMAX/PRAD 
GAIII0IE=4*PI*UMAII0IE/PRADI0IE 
GDB= 1 0 . *al og 10 ( G All ) 

GDBIQIE=10 . ealoglO(GAIIIOIE) 

HRITEC6,*) ’GAII’ .IDIVT, ’  =  '.GAII,’  II  DB  *  ’ ,GDB 
WRITE(6,*) ’ GAIIIOIE  =  ’.GAIIIOIE,'  II  DB  =  ' , GDBIOIE 
WRITE(7,*) ’IDIVT( ' .IDIVT, * )*’ .IDIVT, ’ ;  IDIVP( * , ID I VP, ’ )*’ .IDIVP 
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WRITE (7 ,2003) ’GAIK  \NDIVT, ’)=’ .GAIN, *• GAIMDB( ’ .NDIVT, ■)= 1 ,GDB 
WRITE(7,2003)’GAIMN0NE(’ , NDIVT, ’)  =  ’.GAIHMOME, 

A ’ ; GAINNONEDB( ’ , NDIVT , ’ )  =  ’  .GDBMONE 

2003  F0RMAT(/ (A , 13, A,F12 . 5.A, I3.A.F10. 5) ) 

C500  continue 
STOP 

999  end 

0*******************************************************,*,*£gQ  OF  MAIN  PROGRAM 
C*** ******** ****** **************** *******************5uBR0UTIHE  PLAIE 
C  REFERENCE :  TECHNICAL  REPORT  TR-80-1  ; (pages  57-64) 

SUBROUTINE  PLANE(M1 , M2 , NP , NT , RH , ZH , XT , AT , THR , R) 

C 

C  PLANE  WAVE  EXCITATION  VECTOR  IN  THE  FAR  FIELD.  FROM  MAUTZ  AND  HARRINGTON. 

C  MODIFIED  TO  DO  ONLY  ONE  ANGLE  AND  FREQUENCY  PER  CALL. 

C  **  NEW  BESSEL  FUNCTION  ROUTINE  FROM  MUM.  RECIPES  ** 

C 

COMPLEX  R( 1600) , U , U1 , UA , UB , FA ( 1 500 ) , FB ( 1 500) , F2A , F2B , FI A , FIB 
COMPLEX  U2.U3.U4.US 

DIMENSION  RH(400) ,ZH(400) ,XT(4) ,AT(4) ,R2(4) ,Z2(4) 

MP=MP-1 

MT=MP-1 

M=HT+MP 

M2=2*M 

CC=COS(THR) 

SS=SIM(THR) 

U*(0.,1.) 

U1=3.141593*U**M1 

M3=M1+1 

M4=M2+3 

IF(Ml.EQ.O)  M3=2 

M5=Ml+2 

M6=M2+2 

DO  12  IP=1 ,MP 

K2=IP 

I=IP+1 

DR= . 5* ( RH ( I ) -RH ( IP ) ) 

DZ=.5*(ZH(I)-ZH(IP) ) 

D1=SQRT(DR*DR+DZ*DZ) 

R1=.25*(RH(I)+RH(IP)) 

IF(Rl.Eq.O.)  Rl=l . 

Zl= . 5*(ZH(I)+ZH(IP) ) 

DR=.5*DR 
D2=DR/R1 
DO  13  L= 1 , NT 
R2(L)=R1+DR*XT(L) 

Z2(L)=Z1+DZ*XT(L) 

13  CONTINUE 
D3=DR*CC 
D4=-DZ*SS 
D5*D1*CC 
DO  23  M=M3,M4 
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FA(M)=0 . 

FB(M)=0 . 

23  COITIIUE 
DO  IS  L=1,*T 
X=SS*R2(L)*2 . 

ARG=Z2(L)*CC 

UA= AT ( L ) *CMPLX (COS ( ARG ) , S II ( ARG ) ) 

UB=XT(L)*UA 

C  THIS  LIME  REPUCES  HARRIIGTOI’S  BLOCK 
DO  25  M=M3,M4 
BES=BESSJ(M-2,X) 

FA(M)=BES*UA+FA(M) 

FB(H)=BES*UB+FB(M) 

25  COITIIUE 
15  COITIIUE 

IF(Ml.IE.O)  GO  TO  26 
FA(1)=-FA(3) 

FB(1)=-FB(3) 

26  UA=U1 

DO  27  M=M5,H6 

M7=M-1 

M8=M+1 

F2A=UA*(FA(M8)+FA(M7) ) 

F2B=UA*(FB(H8)+FB(H7) ) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7) ) 

F1B=UB*(FB(H8)-FB(H7) ) 

U4=D4*UA 

U2=D3*F1A+U4*FA(H) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+I 

K5=K2+I 

R ( K2+KT ) =-D5* ( F2 A+D2*F2B ) 

R(K5+HT)=D1*(F1A+D2*F1B) 

IF(IP.EQ.l)  GO  TO  21 
R(K1)=R(K1)+U2-U3 
R(K4)=R(K4)+U4-U6 
IF(IP.EQ.MP)  GO  TO  22 

21  R(K2)=U2+U3 
R(KS)=U4+U6 

22  K2=K2+I2 
UA=UB 

27  COITIIUE 
12  COITIIUE 

RETURI 

BID 

O****************************************SUBR0UTIIE  AITFF 
SUBROUTINE  AITFF (THETA , PHI , THETAS , PHIS , ARAD , ETF , EPF ) 


C  COMPUTES  CLOSED  FORM  AITENNA  PATTER*  II  THE  FAR  FIELD. 

C  ETF  AND  EPF  ARE  RETURNED.  ANGLES  ARE  IN  RADIANS. 

C  PRESENT  ENTRY  IS  FOR: 

C 

C  UNIFORM  ILLUMINATION  SCANNED  TO  A  DIRECTION  (THETAS, PHIS) 

C 

COMPLEX  ETF , EPF , JK , SCL , EB 
JK=(0. 0,6. 283185307) 

BK=6 . 283185307 
PI=BK/2 . 

C  SET  RADIATIO*  PATTERN  TO  ZERO  IN  THE  REAR  HEMISPHERE 
IF(COS(THETA) .LT.O . )  THEN 
EB=(0 . ,0.) 

ELSE 

BB=SIH(THETA)*SIN(PHI)-SIN(THETAS)*SII(PHIS) 

AA=SIN (THETA ) *COS(PHI ) -SIN (THETAS ) *COS (PHIS ) 
ARG=BK*ARAD*SQRT(AA**2+BB**2) 

SCL=-(0. ,1. )*2.*PI**2 
IF(ABS(ARG) .LT. 1 .E-5)THE* 

EB=SCL/2. 

ELSE 

EB=SCL* (BESS J 1 ( ARG ) / ARG ) 

ENDIF 

ENDIF 

C  EB  IS  THE  X  COMPONENT ;  GET  ETHETA  (ETF)  AID  EPHI  (EPF) 

C  ASSUME  AI  ELEMENT  FACTOR,  ELMT 
ELMT=SQRT (ABS(COS( THETA ) ) ) 

ETF=EB*COS (THETA ) *COS (PHI ) *ELMT 
EPF=-EB*SII (PHI ) *ELMT 
RETUR* 

END 

C#*#*#*****************************************BESSEL  FUICTIOIS 
FUNCTIO*  BESSJ(N.X) 

C  RETURNS  THE  BESSEL  FUICTIOI  B  OF  ORDER  I  (>1)  AID  REAL 
C  ARGUMENT  X. 

PARAMETER  (IACC=40 , BIGI0=1 . E10 , BIGNI=1 . E- 10) 

IF(N.EQ.O)  THEI 
C  IF  1=0  CALL  BESSJO 

BESSJ=BESSJO(X) 

ELSE  IF(I.Eq.l)  THEI 
C  IF  *=1  CALL  BESSJ1 

BESSJ=BESSJ1(X) 

ELSE 

C  IF  I>1  USE  RECURS I 01 
BESSJ-O. 

IF(X.IE.O. )  THEI 
T0X=2./X 

IF(X.GT.FL0AT(I))  THEI 
BJMsBESSJO(X) 

BJ=BESSJ1 (X) 

DO  11  J*1,I-1 


BJP=J*TOX*BJ-BJH 

BJM=BJ 

BJ=BJP 

11  COITIIUE 
BESSJ=BJ 
ELSE 

M=2* ( (B+IIT (SQRT (FLOAT( IACC*I ) ) ) )/2) 

BESSJ=0 . 

JSUM=0. 

SUM=0 . 

BJP=0. 

BJ=1 . 

DO  12  J=M, 1 , -1 

BJM=J*TOX*BJ-BJP 

BJP=BJ 

BJ=BJH 

IF(ABS(BJ) .GT.BIGIO)  THEM 

BJ=BJ*BIGII 

BJP=BJP*BIGII 

BESS J=BESS J*BIGII 

SUM=SUM*BIGII 

EIDIF 

IF(JSUM.IE.O)  SUM=SUH+BJ 
JSUM=1-JSUM 
IF(J.EQ.I)  BESSJ=BJP 

12  COITIIUE 
SUM=2.*SUM-BJ 
BESSJ=BESSJ/SUM 
EIDIF 

EIDIF 

EIDIF 

RETURI 

EID 

FUICTIM  BESSJO(X) 

C 

C  BESSEL  FUICTIOI  OF  0  ORDER,  REAL  ARGUHEIT  X 
C  (SEE  'IUMERICAL  RECIPES',  P.172) 

C 

REAL*8  Y,P1 ,P2,P3,P4,P5,01 ,Q2,Q3,Q4,Q5,R1 ,R2,R3,R4,R5,R6, 

*  S 1 , S2 , S3 , S4 , S6 , S6 

DATA  PI ,P2 , P3 ,P4 , P5/ 1 . DO , - . 109862827D-2 , . 27345 10407D-4 , 

*  -.2073370839D-S, .20938872110-6/ 

DATA  qi ,q2 ,Q3 ,Q4 ,QS/- . 186249999SD- 1 , . 143048876SD-3 , 

*  -.6911147661D-5, .782109S161D-6, -.9349451620-7/ 

DATA  R1,R2,R3,R4,R5, R6/57568490574 . DO , - 13362590354 . DO , 

*  651619640. 7D0, -11214424. 18D0, 77392. 33017D0, -184. 9062456D0/ 
DATA  S 1 , S2 , S3 , S4 , S5 , S6/575684904 11.00,1 029632985 . DO , 

*  9494680 . 7 18D0 , 59272 . 64853D0 , 267 . 8532712D0 , 1 . DO/ 

IF(X.Eq.O.)  THE! 

BESS JO* 1 . 

ELSE  IF(ABS(X).LT.8.)  THEM 
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Y=X**2 

BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6) ) ) ) )/ 

*  (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ) )) 

ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX-. 785398164 

BESSJO=SqRT(. 636619772/ AX)*(C0S(XX)*(P1+Y*(P2+Y*(P3+ 

*  Y*(P4+Y*PS) ) ) )-Z*SII(XX)*(Ql+Y*(Q2+Y*(Q3+ 

*  Y*(Q4+Y*Q5))))) 

EMDIF 

RETURH 

EID 

FUICTIOI  BESSJl(X) 

C 

C  BESSEL  FUICTIOI  B  OF  ORDER  1,  REAL  ARGUMEIT  X 
C  (SEE  'IUMERICAL  RECIPES',  P.173) 

C 

REAL*8  Y,P1 ,P2,P3,P4,P5,Q1 ,Q2,Q3,Q4,Q5,R1 ,R2,R3,R4,R5,R6, 

*  S 1 , S2 , S3 , S4 , S5 , S6 

DATA  PI ,P2,P3 ,P4,P5/1 .DO, . 183105D-2 ,- . 3516396496D-4 , 

*  . 2467520174D-5 , - . 20337019D-6/ 

DATA  Q 1 , Q2 , q3 , Q4 , Q5/ . 0468749999SD0 , - . 2002690873D-3 , 

*  . 8449199096D-5 , - . 99228987D-6 , . 10S787412D-6/ 

DATA  R1,R2,R3,R4,R5, R6/723626 14232 . DO , -7896059235 . DO , 

*  242396853 . IDO ,-2972611 . 439D0 , 15704. 4826D0, -30- 16036608D0/ 
DATA  S 1 , S2 , S3 , S4 , S5 , S6/144725228442 . DO , 2300535 178 . DO , 

*  18583304 . 74D0 , 99447 . 43394D0 , 376 . 9991397D0 , 1 . DO/ 
IF(X.Eq.O.)  THEI 

BESS J 1=0. 

ELSE  IF(ABS(X) .LT.8. )  THEI 
Y=X**2 

BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/ 

*  (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ) ) ) 

ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX-2. 356 194491 

BESSJl=SqRT( ,636619772/AX)*(C0S(XX)*(P1+Y*(P2+Y*(P3+ 

*  Y* ( P4+Y*P5 ) ) ) ) -Z*SII (XX )*(qi+Y*(q2+Y*(q3+ 

*  Y*(q4+Y*q5)))))*SIGI(l.,X) 

EID  IF 

RETURI 

EID 


C***********; 
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COMMEITS  :  THIS  SUBROUTIIE  HILL  RIGOROUSLY  CALCULATE  THE  ELECTRIC  FIELD 
FOR  A  CIRCULAR  APERTURE.  THE  FIELD  IS  CALCULATED  AT  THE  COORDIIATES 
SPECIFIED  BY  RH( . )  A HD  ZH( . )  THE  APERTURE  IS  LOCATED  AT  Z  =  0, 

AHD  SCAIIED  TO  A  DIRECTIOH  (THS.PHS).  THE  SUBROUTIHE  OGIVE  IS  THE 
SOURCE  OF  THE  GEOMETRIC  DATA  REQUIRED  BY  CIRC  TO  PERFORM  COMPUTATIOIS . 
ALL  PHYSICAL  DIMEISIOHS  ARE  IORMALIZED  TO  WAVELEHGTH . 


C************* ******************************* ********************************* 


c 

+ant 
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_ 

_ 
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c 
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.  ..)*TAPER(rp) 

c 
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o 
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>% 
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* 

* 

* 
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c***************************************************************************** 

SUBROUTIHE  CIRCRTP ( CHPHI , XP , AP , ARAD .THS.PHS, 

*  PHI , RZ , Z , CIRCR , CIRCT , CIRCP) 

IHTEGER  CHPHI 
DIMEHSIOH  XP(100),AP(100) 

REAL  K 

COMPLEX  J, JK.Gl ,G2,X1 ,COH,CC,SUMX,SUMY,SUMZ 
COMPLEX  CIRCR. CIRCT. CIRCP 
PI=3. 141592654  • 

ETA=120.*PI 

C  SET  FIELD  COMPOHEHTS  TO  ZERO  II  THE  REAR  HEMISPHERE 
CIRCR=(0. ,0.) 

CIRCT=(0. ,0. ) 

CIRCP=(0. ,0.) 

C  Z  OF  THE  OBSERVATIOI  POIHT  IS  >=  0  II  THE  FORWARD  HEMISPHERE 
IF(Z.GE.O.)  THE! 

SUMX=(0. 0,0.0) 

SUMY=(0. 0,0.0) 

SUMZ=(0. 0,0.0) 

K=6. 283185307 
J  *  (0.0, 1.0) 

JK=  (0.0,6.283185307) 

C0I=-J/(4.*PI) 

C  OMIT  THE  FACTOR  OF  ETA  SIICE  IT  IS  HOT  IICLUDED  II  SUBROUTIIE  ZMAT. 

C  ARAD**2  IS  FOR  SCALIIG  OF  GAUSSIAI  IITEGRATIOI .  READJUST  SCALIIG 
C  TO  AGREE  WITH  OLD  SUBROUTIIE  SPHERE  (ARBITRARY) 

CC=  COI 

C  OUTER  LOOP:  IITEGRATE  OVER  X,  -ARAD  <  X  <  ARAD. 

C  IIIER  LOOP:  IITEGRATE  OVER  Y,  -ARAD  <  Y  <  ARAD. 

C  EVALUATE  AITEIIA  FIELD  AT  RZ,Z: 

XO»RZ*COS(PHI) 

YO»RZ*SII(PHI) 
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STS=SI*(-THS) 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 

C  IITEGRATE  II  X,Y  COORDIIATES 
DO  50  M=1 ,CIPHI 
X=ARAD*XP(M) 

DO  60  1=1 ,CIPHI 
Y=ARAD*XP(I) 

RP=SQRT(X**2+Y**2)  « 

C  SKIP  IHTEGRATIOI  POUTS  BEYOKD  THE  AHTEHIA  RADIUS 
IF(RP.LE.ARAD)  THEM 
PHIPRIME=ATAI2(Y , X+ 1 . E- 10 ) 

R=SqRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIH((PHI-PHIPRIME)/2)**2)) 

G1=(((K*R)**2)-1-(JK*R))/R**3 

G2=(3+(3*JK*R)-(K*R)**2)/R**5 

X1=JK*(RP*C0S(PHIPRIME-PHS)*STS-R) 

Y1=X0-X 

X2=Y1**2 

Y2=Y0-Y 

SUKX=SUMX+(CC^AP(I!)^AP(H)*(G1+(X2*G2))*CEXP(X1))*TAPER(RP) 

SUMY=SUMY+(CC*AP(M)*AP(I) *Y1*Y2*G2*CEXP (XI ) ) *TAPER(RP) 

SUMZ=SUMZ+ ( CC*AP (M ) * AP ( I ) *Y1 *Z*G2*CEXP (XI)) *TAPER ( RP ) 

EIDIF 

60  COITIIUE 
50  COITIIUE 

C  ASSUME  Af  ELEMEKT  FACTOR,  ELMT 
ELMT=SqRT(ABS(C) ) 

CIRCT= ( C*COS (PHI ) *SUMX+C*SII (PHI) ♦SUMY-S*SUMZ ) *ELMT 

CIRCR=(S*COS(PHI)*SUMX+S*SII(PHI)*SUMY+C*SUMZ)*ELMT 

CIRCP= ( -SII (PHI ) ♦SUMX+COS (PHI ) *SUMY) ♦ELMT 

EIDIF 

RETURI 

EID 

^♦♦♦♦••••♦♦•♦•♦♦♦♦♦♦♦♦♦♦********»*********FUICTIOI  TAPER 
FUICTIOI  TAPER(RHO) 

C  SPECIFY  AITEIIA  ILLUMIIATIOI  FUICTIOI.  REAL  FUICTIOI  OF 
C  RHO  OILY 

TAPER=1 . 

RETURI 

EID 


f- 
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